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ABSTRACT 


The  regenerative  method  for  stochastic  system  simulation 
allows  data  collection  each  time  the  stochastic  process  enters 
a specific  single  state , r , called  the  regeneration  state. 

The  generated  observations  have  the  desireable  property  of  being 
independent  and  identically  distributed.  Relative  to  a fixed 
run  length,  however,  the  mean  time  between  entries  into  r may 
be  excessively  long  for  complicated  stochastic  systems,  thus  pro- 
viding few  observations  and  poor  variance  estimates. 

The  almost  regenerative  method  is  an  extension  of  the  regener- 
ative method  designed  to  alleviate  this  problem  for  complicated 
stochastic  systems  (such  as  a network  of  queues) . The  almost 
regenerative  method  allows  data  collection  each  time  the  stochastic 
process  enters  a set  of  states.  Simulations  of  simple  queueing 
networks  show  that  the  almost  regenerative  method  can  provide  an 
order  of  magnitude  improvement  over  the  regenerative  method  in 
terms  of  the  mean-square-error  of  the  estimator  of  total  delay 
in  queue,  and  this  relative  improvement  increases  with  system 
complexity.  In  addition,  variance  estimates  based  on  observations 
generated  by  the  almost  regenerative  method  are  shown  emperically 
to  be  superior  to  variance  estimates  based  on  observations  generated 
by  the  regenerative  method. 

The  almost  regenerative  method  is  also  shown  emperically  to 
consistently  provide  mere  efficient  estimators,  in  terms  of  mean- 
square-error,  than  the  frequently  used  fixed  time  increment  method 
for  event-oriented  simulations. 

The  computer  program  modifications  necessary  to  implement  the 
almost  regenerative  method  are  no  more  complicated  than  implementing 
the  regenerative  method  or  the  fixed  time  increment  method.  Also, 
CPU  requirements  are  virtually  the  same  for  each  of  these  methods. 


TABLE  OF  CONTENTS 


l 


k 

I 


I 

i 


i 

i 


i 


Page 


Chapter  1:  Introduction  and  Sumuwry  of  Reaulta  ...  1 

1.1  Introduction  1 

1.2  bummary  of  Results 3 

Chapter  2:  The  Almost  Regenerative  Method  . 6 

2.1  Regenerative  Processes  and  the  Regenerative  Method  (RM).  ...  6 

2.2  The  Almost  Regenerative  Method  (aRM)  . 11 

2.3  The  Asymptotic  Distribution  of  the  Almost  Regenerative 

Estimator 19 

2.4  The  Ratio  Estimators 23 

Chapter  3:  Application  of  the  Almost  Regenerative  Method  to 

a Tandem  Queueing  Network  28 

3.1  The  Almost  Regenerative  State  Space  Approach  28 

3.2  Simulation  of  the  M/M/1  Tandem  Network  . 32 

3.3  Tests  of  the  Renewal  Hypothesis 40 

3.4  Distribution  of  the  Times  Between  Transitions 42 


Chapter  4:  Using  the  Almost  Regenerative  Method  to  Obtain 

Order  of  Magnitude  Improvements  in  MSE  Efficiency  ....  49 


4.1  The  M/M/c  Tandem  Queue  Network  49 

4.2  Simulation  of  the  M/M/2  Tandem  Network  Using  ARM(A(n.  ^nj)]*  • • 51 

4.3  Comments  on  the  Implied  Power  of  the  Almost  Regenerative 

Method  59 


1 


I 

i 

I 

I 


Chapter  5:  The  Almost  Regenerative  Method  in  Perspective  61 

Appendix  I:  Testing  for  Renewal  Processes 64 

References 72 

List  of  Tables: 

Table  3.1:  M/M/1  Tandem  Queue:  Comparison  of  MSE  Estimators*  * 36 

Table  3.2:  M/M/1  Tandem  Queue:  Comparison  of  Variance 

Estimators  .....  36 

Table  3.3:  H^/H^/l  Tandem  Queue:  Comparison  of  MSE  Estimators.  39 

Table  3.4:  H^/H^/l  Tandem  Queue:  Comparison  of  Variance 

Estimators 39 

Table  3.5:  Summary  of  Renewal  Test  Statistics  for  the  Tandem 

Queue  Model  ...... 44 

Table  3.6:  Selected  Values  of  CV  tl|2n}  , CVZ{T|2n+l}  , and 

CV2{T} 48 

Table  4.1:  M/M/2  Tandem  Queue:  Comparison  of  MSE  Estimators*  • 54 

Table  4.2:  M/M/2  Tandem  Network  Comparison  of  MSE  Efficiencies*  55 

Table  4.3:  M/M/2  Tandem  Queue:  Comparison  of  Variance 

Estimators 56 


TABLE  OF  CONTENTS 


List  of  Tables  cont. 

Table  4.4:  M/M/2  Tandem  Queue:  Renewal  Test  Statistics  for 


ARM(A(n1,n2)]  57 

Table  4.5:  Two  M/M/c  Queues  in  Tandem:  Comparison  of  E{T} 

for  the  ARM  and  the  KM 60 


List  of  Figures: 


Figure  3.1:  Autocorrelation  Function  for  25  Lags:  M/M/1 

Tandcn  ARM(A(n)] 41 

Figure  3.2:  Autocorrelation  Function  for  25  Lags:  M/M/1 

Tandem  FTM 41 

Figure  3.3:  Autocorrelation  Function  for  25  Lags:  H./H./l 

Tandem  ARM[A(n)] 43 

Figure  3.4:  Autocorrelation  Function  for  25  Lags:  H./H./l 

Tandem  FTM . . . . 43 


CHAPTER  1 

| INTRODUCTION  AND  SUMMARY  OF  RESULTS 

1.1  Introduction 

The  central  problem  in  simulations  of  stochastic  processes  is 
concerned  with  the  method  used  to  collect  observations  of  the  sto- 
chastic process  und.  r consideration.  Consider  a stationary  stochastic 
process  {X(t),  t > 0}  , and  let  f be  a real-valued,  integrable 
function  defined  on  X(t)  . For  a fixed  simulation  run  length  £ , 
an  estimate  of  6 - E{f(X)}  is 

6(£)  =7  / f [X(u)]du  , 

*\) 

| 

j since  it  can  be  shown,  under  fairly  general  conditions,  that  0(0  -*■  0 

l with  probability  one  (w.p.l)  as  £•*■“.  Thus,  observations  nay  be 

I 

i taken  at  arbitrary  intervals  to  form  an  estimate  of  0 . However,  the 

i 

j 

l conditions  under  which  observations  are  collected  is  important  when 

estimating  the  variance  of  6(0  (which  is  used,  for  example,  in  con- 
fidence interval  estimation  and  hypothesis  testing),  since  the  obser- 
vations may  be  serially  correlated.  All  of  the  methods  considered  in 

i 

this  dissertation  partition  the  time  series  (X(t),  t c [0,£]}  in 

i 

i 

different  ways  in  an  effort  to  reduce,  or  eliminate,  the  dependence 
between  observations.  The  difference  between  the  methods  depends  on 
the  use  made  of  knowleege  about  the  process  X(t)  . 

The  classical  j'ixed  tine  increment  method  (FTM)  is  the  degenerate 
case,  since  no  use  is  made  of  any  knowledge  of  X(t)  . Observations 
are  collected  at  fixed  interval  lengths,  say  A , where  nA  * £ and 

I. 
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n is  the  number  of  observations.  Since  these  observations  may  be  de- 
pendent, the  effects  of  autocorrelation  must  be  allowed  for  in  the 

A 

estimate  of  V{0(£)}  » thus  creating  an  estimation  problem.  Two 
approaches  to  this  problem  have  been  the  method  of  batch  means 
[Conway  (1963)]  and  the  application  of  time  series  methods  [Conway, 
Johnson,  and  Maxwell  (1939),  and  Fishman  and  Kiviat  (1967)]. 

First  suggested,  but  not  developed,  by  Cox  and  Smith  (1961),  the 
regenerative  method  (RM)  for  collecting  observations  of  the  time  series 
(X(t),  t e [0,£]}  has  recently  been  systematically  developed  by  a 
number  of  authors  [Kabek  (1968),  Fishman  (1973,  197A  abc) , Crane  and 
Iglehart  (1973,  1974  abc,  1975  ab),  and  Iglehart  (1974,  1975)].  The 
stochastic  process  (X(t),  t > 0}  is  called  regenerative  (a  more 
precise  definition  is  given  in  Chapter  2)  if,  every  time  a certain 
event  occurs,  the  future  of  the  process  from  then  on  is  independent 
of  the  past  and  becomes  a probabilistic  replica  of  the  future  after 
time  zero.  Such  times  are  called  regeneration  points.  The  central 
idea  of  the  RM  is  to  exploit  the  fact,  when  X(t)  is  regenerative, 
that  certain  random  /ariables  between  successive  regeneration  times 
are  independent  and  identically  distributed  (iid),  thus  circumventing 
the  autocorrelation  problem  in  estimates  of  V{0(£)}  . 

As  a practical  issue,  the  stochastic  process  X(t)  may  not  have 
the  regenerative  property  (e.g.,  a stable  multiple  channel  queue  with 
interarrival  and  service  time  distributions  such  that  the  system  never 
empties  [Whitt  (1972)]),  or  if  it  does,  the  expected  time  between  re- 
generation times  may  be  excessively  large  (e.g.,  a network  of  queues). 
In  either  of  these  cases  it  may  be  possible  to  define  an  almost  re- 
generative method  (ARM)  which  retains  the  dcsireable  properties  of 
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the  RM,  but  which  allows  more  frequent  data  collection  than  might 

be  possible  if  the  RM  were  applied  directly. 

To  see  that  sample  size  is  important  when  estimating  V{ 0 } 

consider  a sequence  of  iid  random  variables  {Z^}J  • each  of  which 

2 

has  a Normal  distribution  with  mean  y and  variance  a ; i.e., 

0 a I** 

Z ^ — N(y ,o  ; for  all  i ■ 1,  . . . , n . Then  0 ■ — ^ 2^  is  an 

A 2 

unbiased  estimator  of  y , and  V{6}  ■ o /n  . The  usual  unbiased 

estimator  of  o is  s = — ^ (Z^  - 0)  . For  samples  from  a 

Normal  distribution,  it  is  veil-known  that  the  sample  function 
o 2 2 

(n  - l)s‘  /a  has  the  X distribution  with  (n  - 1)  degrees  of 

freedom.  Thus,  for  any  probability  1 - a , one  can  select,  two  values 
2 

of  X , depending  on  n and  a , such  that 


(1.1) 


a . 


As  a result,  this  confidence  interval  can  be  made  as  small  as  desired 
simply  by  increasing  n . 

The  purpose  of  this  dissertation  has  been  to  investigate  and 
develop  the  properties  of  an  almost  regenerative  method  suitable 
for  application  to  variance  estimation  in  complicated,  practical 
queueing  systems  which  do  not  have  frequently  occuring  regeneration 
points. 


1.2  Summary  of  Results 

The  regenerative  method  (KM)  and  the  almost  regenerative  method 
(ARM)  are  defined  and  discussed  in  Chapter  2.  The  ARM  estimators  are 
found  to  have  the  same  asymptotic  properties  as  the  RM  estimators  when  the 
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stochastic  process  under  consideration  is  regenerative. 


When  applying  the  RM  to  a stochastic  process  {X(t),  t > 0)  with 


state  space  E , the  main  idea  is  to  select  a single  state  in  E such 


that  the  X(t)  process  regenerates  each  time  it  enters  this  state.  The 


ARM  permits  one  to  define  a set  of  states , say  A , of  E such  that  the 


X(t)  process  "almost  regenerates,"  with  negligible  dependence  on  the 


past  history  of  the  process,  each  time  X(t)  enters  A . As  a result, 


the  ARM  will  produce  more  observations  than  the  RM  for  a fixed  run 


length.  A good  choice  of  A depends  on  the  model  parameters  being 


estimated. 


The  observations  generated  by  the  ARM  for  the  queueing  networks 


simulated  in  Chapters  3 and  4 are  found  to  act  as  if  they  were  iidt 


thus  allowing  direct  use  of  the  regenerative  estimators. 


A natural  way  to  compare  the  relative  efficiencies  of  the  ARM, 


RM,  and  FTM  estimators  is  to  compare  their  mean-square-errors,  where 


the  mean-square-error  of  the  estimator  6 of  6 is  defined  by 


MSE  {6}  ■ E{  (6  - 0)  } . The  ARM  estimators  for  the  queueing  network 


models  of  Chapters  3 and  4 were  found  emperically  to  be  up  to  60  percent 


more  efficient  than  either  the  RM  or  FTM  estimators.  This  result 


is  attributable  to  the  facts  that  the  ARM  can  generate  observations 


which  act  iid,  unlike  the  FTM,  and  the  ARM  can  generate  more  ob- 


servations than  the  RM  for  the  same  run  lengths.  Examples  are  given 


where  one  might  expect  the  improvement  in  efficiency  of  the  ARM 


estimators  to  be  up  to  an  order  of  magnitude,  or  more,  more  efficient 


than  comparable  RM  or  FTM  estimators. 


When  estimating  the  variance  of  the  delay  in  queue  estimator, 


the  ARM  variance  estimators  were  found  emperically  to  be  superior 


_2_ 
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to  the  RM  variance  estimators;  i.e.,  the  ARM  variance  estimators 
have  a smaller  bias  and  a smaller  estimated  standard  deviation  than 
the  RM  variance  estimators. 

Additionally,  comparisons  of  the  simple  ratio  estimator  with  the 
jackknife  ratio  estimator  and  comparisons  of  the  simple  variance 
estimator  with  the  jackknife  variance  estimator  show  emperically 
that  the  jackknife  estimators  are  preferable  for  ARM  generated  obser- 
vations. 
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CHAPTER  2 

THE  ALMOST  REGENERATIVE  METHOD 


The  regenerative  method  for  simulating  stochastic  systems  Is 
reviewed  In  this  chapter,  and  an  extension  of  this  method,  the  almost 
regenerative  method,  Is  presented  which  is  applicable  to  simulations 
of  large,  complicated  stochastic  systems.  In  Chapters  3 and  A, 
examples  of  queueing  networks  are  investigated  which  show  emperically 
that  the  almost  regenerative  estimators  can  be  significantly  more 
efficient  than  the  regenerative  estimators. 

2.1  Regenerative  Processes  and  the  Regenerative  Method  (Ril) 

A sequence  Y^^  •••  of  non-negative,  independent  and  identic- 
ally distributed  (iid)  random  variables  is  called  an  ordinary 
renewal  process.  To  avoid  trivialities,  assume  that  P{Yn  > 0}  > 0 
for  all  n . The  physical  principle  underlying  renewal  theory  is 

that  of  events  occuring  at  arbitrary  epochs,  where  the  intervals 

between  these  events  are  the  Y . Define  S„  ■ 0 and  S - Y + S , 

n 0 n n n-1 

n ■ 1,2,  ....  Then  Sr  is  the  epoch  of  the  nth  renewal  alter 

epoch  0 . Wien  Y^  has  a different  distribution  than  the  other 

Yq,  n > 1 (Y^,Y2»  ...  are  still  independent),  the  sequence  Y^,Yj  ... 

is  called  a general  (delayed)  renewal  process. 

A stochastic  process  {X(t),  t > 0}  , with  state  space  E , 

is  said  tp  be  a regenerative  process  if  there  exists  a sequence  of 

random  variables  {Y  , with  S_  - 0 and  S ■ Y + S , , 

n 1 0 n n n-1  * 

n ■ 1,2,  ...  , such  that 

OS 

(Rl)  {Y^}^  is  an  ordinary  renewal  process;  and, 

(R2)  for  any  n ■ 0,1,  ...  , and  m ■ 1,2,  ...  , 
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P{X(Sn  + t£)  c Ej,  i - 1,2 m | X(«),  ■ < Sn)  ■ 

p{x(s0  + tA)  e EA,  i - 1,2,  ....  n)  , 

where  Ei  C E and  r > 0 , for  i " 1,2,  ...,  m . The  epochs 
S^.Sj,  ...  are  called  regeneration  points  (epochs),  and  the 
Y^.Y^,...  are  the  lengths  of  the  regenerative  cycles.  Thus,  a 
regenerative  process  is  one  which  "starts  over"  at  each  regeneration 
point  in  such  a way  that  the  future  of  the  process  at  a regeneration 
point  is  independent  of  the  past  history  of  the  process,  and,  further- 
more, the  future  at  a regeneration  point  is  a probabilistic  replica 

00 

of  the  future  after  epoch  . When  the  {Yn}^  , form  a general 

renewal  process,  { X ( t ) , t > 0}  is  said  to  be  a general  (delated) 

regenerative  grocers.  In  this  case,  Sq  is  replaced  by  in 

(R2).  Unless  specified  to  the  contrary,  all  regenerative  processes 

considered  in  this  chapter  will  be  ordinary  regenerative  processes. 

So  far  nothing  has  been  said  about  the  state  of  X(t)  at  the 

00 

regeneration  points  { }q  . It  is  possible  that  a state  change 
in  the  X(t)  process  may  not  accompany  the  regeneration  points. 

For  example,  although  the  GI/G/1  queue  is  regenerative  ..  those 
epochs  when  an  arrival  finds  no  customers  in  the  system,  the  number 
of  customers  in  queue  does  not  change  (it  remains  zero).  It  is 
also  possible  that  state  changes  in  X(t)  at  regeneration  points 
nay  not  be  union Consider  the  number  in  system  process  tor  the 
GI/G/1  queue  witli  batciies  of  random  size.  At  the  epoch  an  arriving 
batch  finds  the  system  empty,  a regeneration  point  occurs  and  the 
number  in  system  becomes  the  batch  size,  which  is  a random  variable. 
Finally,  for  a given  X(t)  process,  there  may  be  n any  choices  of 
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regeneration  points.  For  example,  the  number  in  system  process  for 

an  M/M/1  queue  has  regeneration  points  when  an  arrival  finds  n 

customers  in  the  system,  for  fixed  n , it  also  has  regeneration 

points  when  a departure  leaves  n customers  behind  in  the  system, 

for  fixed  n . This  property  is  true  for  all  n ■ 0,1,  .... 

To  avoid  later  confusion  on  the  above  issues,  in  this  dissertation 

X(t)  will  represent  a stochastic  process  such  that,  if  X(Sq)  - u 

and  XCS^)  - v (i.e.,  the  states  of  X(t)  just  before  and  just 

after  the  nth  r/egeneration  point)  for  some  fixed  n and  u,v  c E , 

then  X(S  ) r-  u and  X(S  ) ■ v for  all  n - 0,1,  ....  The  state 
n n 

pair  (u,v)  will  be  referred  to  as  the  regenerative  state  transition 
for  the  profess  X(t)  , and  X(t)  will  be  referred  to  as  a (usv) 
regenerative  process  whenever  necessary  to  distinguish  the  states 
u and  v . This  terminology  is  convenient  and  compact;  it  will 
become  clearer  when  "almost  regenerative  processes"  are  introduced 
in  which  the  state  transitions  of  interest  are  allowed  between  sets 
of  states.  Further,  this  terminology  will  not  restrict  the  generality 
of  our  results. 

Now  consider  a regenerative  process  (X(t),  t > 0}  which 

converges  in  distribution  to  a random  variable  X , denoted  by 

, . V 

X(t)  X , and  suppose  a point  estimate  of  the  constant  0 » E{ f (X) } 
is  desired,  where  f is  a real-valued,  integrable  function  defined 
on  X . For  example,  if  f(x)  ■ x then  0 is  the  mean  of  X . If 
f(x)  - 1 when  x < a and  f(x)  “0  when  x > a , then  0 * P{X  < a)  . 
A cumulative  process  {An>  n - 1,2,  ...  ) may  be  defined  on 


X(t)  by 
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(1.1) 


A 

n 


f [X(u) ]du 


where  the  { } are  the  regeneration  points  of  X(t)  [corresponding 
to  a specific  choice  of  regenerative  state  transition  (u,v)].  Since 
X(t)  is  regenerative,  the  { are  iid.  When  E^1  |f[X(u)]|du  < ® 


[Brown  and  Ross  (1970)], 


t 

(1.2)  lin  f f[X(u)]du/t  - 0 

t xx.  J0 


with  probability  one  (w.p.l).  And  from  Smith  (1955),  Theorem  7, 


(1.3) 


t 

lim  f fiX(u)]du/t 
t-~  J0 


e{a1) 

"e(y \ W*P’1  • 


Relations  (1.2)  and  (1.3)  are  also  true  in  expectation.  From  (1.2) 
and  (1.3)  it  is  immediate  that 


(1.4) 


E{A1J 

" e{yJT  ‘ 


When  X(t)  is  a general  regenerative  process,  this  result  can  be 
restated  as  0 = E{A2}/E{Y2). 

The  form  of  (1.4)  suggests  a straight-forward  method,  called  the 
regenerative  method  (RM) , for  estimating  0 . Suppose  the  regenerative 
process  (X(t),  t > 0}  is  observed  for  n regenerative  cycles  and 
the  pairs  (A^Y^  are  recorded  for  cycle  i , i - 1,2,  ...»  n , 
where  Y^  is  the  ith  cycle  length  and  is  defined  by  (1.1). 


Since  Che  pairs  {(A  ,Y^)}J  are  iid,  it  follows  from  the  Strong  Law 
of  Large  Numbers  (SLLN)  and  (1.4)  that 

(1.5)  lim  7^  - 0 w.p.l  , 

n+®  Y(n) 

where  A(n)  - A^/n  and  Y(n)  - Y^/n  . Thus,  a regenerative 
estimator  of  6 is 

(1.6)  0 (n)  - , n - 1,2 

Y(n) 

The  RM  eliminates  some  of  the  shortcomings  of  the  fixed  time 
increment  method  (FTM) , which  was  described  in  Chapter  1.  The 
observations  generated  by  the  FTM  may  be  correlated,  thus  complicating 
estimates  of  estimator  variance;  there  is  no  autocorrelation  between 
the  RM  observations  {(A^.Y^)}"  which  are  mutually  independent. 

Also,  the  FTM  observations  are  influenced  by  the  initial  state  of  the 
system,  thus  encouraging  one  to  "throw  away"  the  first  few  observations 
to  reduce  estimator  bias  [Conway  (1963)];  sin^e  the  KM  observations 
are  identically  distributed,  data  can  be  collected  from  the  beginning 
of  the  simulation. 

However,  the  RM  does  have  some  drawbacks  of  its  own.  First, 
the  ratio  estimator  (1.6)  of  9 is  biased  because,  in  general,  the 
expected  value  of  a ratio  is  not  the  ratio  of  the  expected  values. 

For  sufficiently  large  sample  sizes,  and  under  the  assumption  that 
the  (A^Y^)  are  iid  bivariate  Normal,  it  is  known  [Tin  (1965)]  that 

E{0(n) } - 0[1  - £ (Cu  - CQ2)]  + 0(n"1) 


(1.7) 
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where 


(1.8) 


'1J 


E{(A1  - EA1)i(Y1  - EY1)j) 
Ei{A1)  EJ{Y1) 


• 1* J ■ 0,1  • . • 


is  the  bivariate  coefficient  of  variation  of  order  (i,j)  , and 

t a|/ 

0(n  ) has  the  property  that  n 0(n  ) < * as  n -*■  * . Note  that 

depends  on  the  choice  of  regenerative  state  transition  (u,v)  . 
Other  ratio  estimators  of  0 which  have  smaller  biases  and  mean- 
square-errors  than  those  of  (1.6)  are  compared  by  Tin  (1965). 

The  KM  has  a more  serious  defect,  however.  The  expected  time 
between  regeneration  points,  E{Y^)  , may  not  be  small  enough  for 
practical  considerations.  Of  course,  it  may  also  be  possible  that 
(X(t),  t > 0}  may  not  be  regenerative,  in  which  case  there  are  no 
states  (u,v)  with  the  desired  property.  In  spite  of  this  problem,  an 
extension  of  the  KM  may  still  be  used  to  collect  data,  from  simulations 
of  the  process,  which  can  be  conveniently  used  to  estimate  6 and  the 
variance  of  the  estimator  of  0 . 

2.2  The  Almost  Regenerative  Method  (ARM) 

The  two  basic  approximation  strategies  one  may  take  when  dealing 
with  stochastic  processes  which  do  not  have  convenient  regeneration 
state  transitions  can  be  classified  as  "state  space  transformation" 
methods  and  "almost  regenerative"  methods. 

The  central  idea  of  the  state  space  transformation  methods 
[Crane  and  Iglehart  (1975b)]  is  to  redefine  the  stochastic  process 
in  such  a way  that  the  new  process  lias  convenient  regeneration  state 
transitions.  The  KM  may  then  be  applied  directly  to  the  new  process, 
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from  which  Inferences  may  be  drawn  about  the  process  of  interest. 

The  danger  with  these  methods  is  that  they  may  produce  estimators 
which  are  not  consistant  in  the  usual  statistical  sense.  For  example, 
consider  a Gl/G/1  queue  with  continuous  interarrival  and  service 
distributions.  By  discretizing  these  two  distributions  to  obtain 
a regenerative  process,  certain  random  variables  may  be  adversely 
affected.  This  can  be  seen  by  considering  the  busy  cycle  length 
process  (the  time  between  successive  arrivals  finding  the  system 
empty).  Since  there  is  a positive  probability  that  the  only  customer 
in  the  system  will  complete  his  service  at  the  same  time  an  arrival 
occurs,  the  definition  of  what  constitutes  a busy  cycle  is  a matter 
of  convention.  If  departures  are  processed  before  arrivals,  the 
new  arrival  begins  a new  busy  cycle,  otherwise  the  new  arrival 
contributes  to  the  continuation  of  the  present  busy  cycle  started 
by  tome  previous  arrival.  Thus,  in  one  case  the  busy  cycles  are 
"short"  and  in  the  other  case  they  are  "long";  neither  may  provide 
a good  approximation  to  the  actual  busy  cycle  length  distribution. 

Unlike  the  state  space  transformation  methods,  the  almost 
regenerative  methods  do  not  change  the  process  of  interest.  Rather, 
they  simply  modify  the  conditions  under  which  observations  are  collected. 
Recall  that  observations  of  the  process,  when  using  the  RM,  are  made 
each  time  the  process  makes  a (u,v)  transition,  and  the  resulting 
observations  are  iid.  This  is  a desireable  property.  As  motivation 
for  the  concept  of  an  almost  regenerative  method,  consider  the  following 
definition  due  to  Smith  (1955),  which  is  a weakening  of  the  definition 
of  a regenerative  process. 

A stochastic  process  (X(t),  t > 0}  with  state  space  E is 


13 


•aid  to  ba  loosely  regenerative  if,  for  n - 0,1,  ...  and  m ■ 1,2,  ... 

thera  exiata  a sequence  of  apocha  (Z  } and  a constant  A = A(E. , ..., 

n i 

E ) > 0 such  that 
m * 

(2.1)  P(X(Z  + t.)  c E.,  i * 1,2 m | X(s),  8 < Z - A} 

nxi  — n 

• P{X(Z  + A + t.)  t E,,  i - 1,2,  ...,  n) 
n 11 

for  E^C  E and  t^  > 0 , i ■ 1,2,  ...,  m . The  difference  between 

a regenerative  process  and  a loosely  regenerative  process  nay  be 

explained  as  follows.  In  a regenerative  process  the  past  history, 

except  in  so  far  as  it  determines  the  most  recent  regeneration 

point,  loses  all  predictive  value  whenever  a regeneration  point 

occurs.  In  a loosely  regenerative  process  the  past  continues  to  have 

some  influence  on  the  behavior  of  the  process  for  a further  amount 

of  time  A after  t ,e  epochs  (Z  } . 

n 

As  an  example  o;'  a loosely  regenerative  process,  consider  a 
GI/G/c  queue  where  service  time  at  each  channel  is  a constant,  say  a , 
and  suppose  the  process  of  interest  is  the  number  of  customers  in 
the  system.  This  process  is  regenerative  at  those  epochs  when  an 
arrival  finds  the  system  empty.  However,  it  is  not  regenerative  at 
those  epochs  when  an  arrival  finds  n(n  > 0)  customers  in  the  system 
since  the  future  behavior  of  the  process  depends  on  the  remaining 
service  time  of  those  customers  in  service,  as  well  as  n . However, 
if  the  history  of  the  process  were  known  for  a tine  units  before 
an  arrival  epoch,  the  remaining  service  times  of  the  customers  in 
service  would  be  known,  and  thus  the  future  probabilistic  behavior 
of  the  process  would  also  be  known.  In  this  case,  the  process  of 
interest  is  loosely  regenerative,  with  A - a , when  an  arrival 
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finds  n customers  in  the  system,  for  fixed  n > 0 . 

This  example  suggests  an  "almost  regenerative"  method  for  simulat- 
ing a queueing  system.  Let  E{T}  be  the  expected  interarrival  time 
of  customers  to  the  above  described  GI/C/c  queue,  and  suppose  that 
the  system  is  stable,  in  the  sense  that  p = u/cE{Tj  » 1 . For 
simulation  purposes,  one  could  take  observations  of  the  process 
each  time  an  arrival  finds  n in  the  systen;  these  epochs  do  not 
constitute  regeneration  points,  but  one  might  expect  the  influence 
of  Lite  remaining  service  times  on  the  future  to  become  negligible 
as  the  time  seperation  between  two  lagged  observations  increases. 

To  see  what  is  meant  by  this  statement,  consider  the  following 
do f ini t ion. 

A stationary  process  (X(t),  t > 0}  with  state  space  E is 
said  to  be  (cf.  Billingsley  (1968)]  if,  for  every  set 

1'  E , there  exists  a monotone  non-increasing,  real-valued  function 
i such  that 

(2.2)  |P{X(t  + B)  t F | X (s) , s < t)  - P{X(t  + 6)  «-  F)  <•  *(fj) 

for  t,fi  > 0 , where  only  functions  A satisfying  lim  i(r!)  - 0 are 

a 

considered.  The  sequence  of  observations  generated  by  an  arrival 
finding  n customers  in  the  system  generally  is  not  stationary,  so 
the  4-mixing  definition  does  not  apply  directly.  However,  the 
essense  of  (2.2)  is  that  separate  observations  of  the  process  are 
asymptotically  independent. 


The  data  collection  method  outlined  above  may  be  generalized. 
Consider  the  stochastic  process  {X ( t ) , t > 0}  with  state  space 
H , and  let  U and  V be  two  disjoint  subsets  of  K (U  and  V 
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need  not  partition  E ).  Also,  let  {Y^}"  be  a sequence  of  non- 
negative random  variables  where  Y^  is  the  time  between  the  (n  - l)st 
and  the  nth  transition  of  X(t)  from  any  state  which  is  an  element 
of  U to  any  state  which  is  an  element  of  V . Define  ■ 0 

and  S'-Y'+s'  , n - 1,2,  ....  The  sequence  { Y ' } _ will  be 

n n n- 1 n l 

called  the  (U,V)  almost  renewal  process  and  (S^q  will  be  called 
the  (U,  V)  almost  regeneration  points  (epochs)  of  the  X(t)  process. 
Note  that  no  conditions  are  placed  on  these  definitions  other  than 
the  transitions  of  interest  be  from  the  set  U to  the  set  V . 

The  following  notational  convention  will  be  observed  throughout 
the  remainder  of  this  chapter:  all  symbols  associated  with  a (U,V) 

representation  of  the  X(t)  process  will  be  so  designated  by  a 
prime  (').  Usually  these  symbols  will  have  a direct  counterpart 
in  a (u,v)  regenerative  process. 

Let  X(t)  5 x and  suppose  a point  estimate  of  the  constant 
0 = E{f(X)}  is  desired,  where  f is  a real-valued,  integrable 
function  defined  on  X . A process  {AM^  may  be  defined  by 


oo 

where  the  {SMq  are  the  (U,V)  almost  regenerative  points  of  X(t) 
corresponding  to  a particular  choice  of  U,V  C E . The  almost 
regenerative  method  (ARM)  for  obtaining  an  estimate  of  0 is 
similar  to  the  regenerative  method  (RM) . Suppose  the  X(t)  process 
is  observed  for  n'  (U,V)  almost  regenerative  cycles,  and  the 
pairs  (A^,Yp  are  recorded  for  cycle  i « 1,2,  ...,  n'  . Then  an 


(2.3) 


S' 

A-  = r 
n Jc , 


f[X(u)]du  , n ■ 1,2 


n-1 


i 
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estimate  of  6 la 

(2.4)  6'(n')  - AW) 

V'(n') 

where  A'  (n‘)  and  Y^n')  are  the  sample  means. 

A 

Conditions  exist  such  that  G'(n')  0 w.p.l  . Although 

n'-K“ 

this  type  of  result  is  felt  to  hold  for  fairly  general  stochastic 

processes,  it  will  be  shown  here  only  when  the  X(t)  process  has 

at  least  one  embedded  (u,v)  regenerative  process;  i.e.,  there 

exist  u , v c E such  that  X(t)  is  (u,v)  regenerative.  Let 

N(n')  be  the  number  of  (u,v)  regeneration  points  during  the  first 

n'  (U,V)  almost  regeneration  points,  and  let  N..  be  the  number  of 

(U,V.)  almost  regeneration  points  during  the  ith  (u,v)  cycle. 

Note  that  the  {N\}  are  iid.  To  avoid  trivialities,  assume  for 

the  remainder  of  this  chapter  that  the  sets  U,V  C E are  chosen  in 

such  a way  that  E{N,  } < °°  and  P{Y’  > 0}  > 0 . 

i n 

The  following  technical  lemma  will  be  needed. 


Lemma  2.1: 

If  X^  is  a random  variable  defined  on  the  ith  (u,v) 
regenerative  cycle  and  E{X^}  < » , then 


(2.5) 


i * "+1  n 

lim  = 0 w.p.i  . 

n  i * * *  v 

n-+® 


Proof : 

Since  the  { X ^ } are  iid,  and  by  t he  SLLN, 


(2.6) 


E{X1> 


n+1 
n + 1 


EfX^  + 11m 
Tt*o° 


II 

Ih 


n+1 


+ lim 


n 

n + 1 


w.p.l.  The  result  then  follows. 

Lemma  2.1  can  easily  be  extended  to  the  following  lemma  which  will 
be  stated  without  proof. 


Lemma  2.2: 

If  X^  ;.s  a random  variable  defined  on  the  ith  (u,v) 
regenerative  cycle,  E{X^}  < « , and  (N(n)}  is  a sequence  of  random 
variables  that  goes  to  « as  n -*•  <*>  , then 


(2.7) 


“ 0 w.p.l  . 


Theorem  2.3: 


If  {X(t),  t > 0}  , with  state  space  E , has  an  embedded 

(u , v)  regenerative  process  such  that  u e U and  v e V for  some 

S 

U,V  CE  , and  if  EJ"n  | f [X  (u)  ] | du  < 00  for  all  n = 1,2,  ...  , 

n-1 

then  the  estimator  (2.4)  is  consistant;  i.e.,  lira  8'(n')  = 6 w.p.l. 


Proof : 


N(n') 

Let  K(n')  = £ N 

1 1 


(2.8)  lim  Y'(n')  = lim 


Then 


ii 

l n 


K(n')  n' 

I V l '■ 

1 J K ( n 1 ) +1 

TL  K(n')  + n'  - K(n’) 





Since  0 < J Y!  < Y„/  , and  since  E{  V . > < » for  all  J 

‘ K(n' )+l  1 ’ N<n  )+1  j 

(Y.  is  the  time  between  the  (j  - l)st  and  the  jth  (u,v)  renewal 
^ n' 

points),  it  follows  from  Lemma  2.2  that  \ Y!/N(n’)  -*■  0 v.p.l. 

K(n')+1  1 

Similarily,  [n'  - K(n')]/N(n’)  -►  0 w.p.l.  Hence,  dividing  both 
numerator  and  denominator  of  (2.8)  by  N(n')  and  applying  the  SLLN 


yields 


(2.9) 


E{Y1) 


lim  Y'(n')  - - , ■ v w.p.l. 

n'~  EtV 


To  obtain  a similar  result  for  A'(n')  , let  the  positive  and 
negative  parts  of  f(x)  be  defined  by  f+(x)  » max  [0,f(x)J  and 
f (x)  = -min  [0, f (x) ] . Then  f(x)  = f+(x)  - f"(x)  for  all  x . 

rs 

Since  E J n |f[X(u)]|du  < 00  for  all  n «*  1,2,  ...  (by  assumption), 
Sn-1  S 

and  this  implies  that  fn  | f £X(u)  ] | du  < « w.p.l  for  all  n = 1,2,  ..., 

Js  ^ 

+ + 

the  same  is  also  true  for  f (x)  and  f (x)  . Define  A',  and  A', 

n'  n 

to  be  the  positive  and  negative  parts  of  A' , . These  terms  are  non- 

n 

negative.  Thus 


(2.10) 


n'  , N(n')+1 

I a;+  < / i 

K(„')+l  s/(nl) 


f [X(u)]du  . 


But  the  right-hand  side  of  (2.10)  has  finite  expectation,  so  application 

n'  + 

of  Lenuna  2.1  yields  T A'  N(n’)  ■*  0 w.p.l  . A similar  result 

K(n' )+l  1 

is  true  for  the  negative  parts  of  A ^ , i = K(n')+1,  ...,  n'  . Hence 
the  same  arguments  that  led  to  (2.9)  yield 


(2.11) 


lim  A' (n' ) 


ECA^ 

E{N±} 


w.p.l  . 
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Then  from  (2.9),  (2.11),  and  (1.4), 

E{A  } 

(2.12)  lim  O'(n')  - - 0 v.p.l, 

n'-~  E{Y1) 

as  was  to  be  shown. 

When  applying  the  ARM  to  a stochastic  process  X(t)  , the 
object  is  to  define  the  sets  U,V  C E such  that  the  serial  correlation 
between  observations  and  the  expected  time  between  successive  (U,V) 
almost  regeneration  points  are  sufficiently  small  to  justify  use  of  the 
regenerative  estimators.  For  example,  to  estimate  the  expected 
delay  in  queue,  d , for  an  M/M/1  model.  Crane  and  Iglehart  (1975b) 
set 

(2.13)  V « {x  : x c (d  - e,  d + c] ) , U - R+  - V 

(where  d was  known)  for  seme  c > 0 . They  found  emperically  that 
values  of  e < d produced  little  autocorrelation  between  successive 
(U,V)  almost  regenerative  cycles.  Interestingly,  they  also  found 
that  as  c was  chosen  to  increase  the  frequency  of  observation, 
the  accuracy  of  the  estimators  actually  improved.  One  would  expect, 
however,  that  this  desireable  property  would  not  persist  as  the 
"trapping  interval"  width  [d  - e,  d + e]  increases  without  bound, 
since  successive  customer  delays  are  correlated.  More  will  be  said 
about  this  issue  in  Chapters  3 and  4. 


2.3  T he  Asymptotic  Distribution  of  the  Almost  Regenerative  Estimator 


It  is  well-known  [cf.  Law  (1974)]  for  the  regenerative  estimator 
, given  by  (1.6),  that  if  0 < E{A^}  < * , 0 < E{Y^}  < 00  , 


0(n)  of  0 


00 
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S 

/n 

| f[X(u)]  | du  < « for  all  n ■ 1,2,  then 

S i 
n- 1 


(3.1) 


Sn  [8(n)  - 8]  V 


/a 


W(0,1) 


n-x» 


where 


V{A,  - 0Y.  ) 

(3.2)  M - — - ez 

EZ{Y  } 


E{A2} 

e2{a1)  e{a1}e{y1)  e2{y1). 


2E{A  Y } E{Y?) 

i— + 


and  W(0,1)  is  the  unit  Normal  distribution  function.  Similar 
results  will  be  shown  in  this  section  for  the  almost  regenera  Live 

A 

estimator  0'(n')  given  by  (2.4),  under  the  condition  that  (X(t), 
t > 0}  is  regenerative. 

The  following  obvious  lemma  will  be  stated  without  proof. 


Lemma  2.4: 

If  the  (X(t),  t > 0}  process,  with  state  space  E , has  at 
least  one  embedded  (u,v)  regenerative  process,  for  some  U,V  c E , 
then  the  expected  number  of  (U,V)  almost  regeneration  points,  , 
during  the  ith  (u,v)  regenerative  cycle  is  at  least  one;  i.e., 


(3.3) 

Theorem  2. 5: 


EfN^  > 1 . 


If  the  {X ( t ) , t > 0}  process  is  as  in  Lemma  2.4,  and  if 

Sn 

0 < E{A^}  < ® , 0 < E{ Y2}  < <*>  , and  hf  |f[X(u)]|du  < °>  for  all 


n-1 


...  . . ....  ■■  - - — 
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n m 1,2,  t • * i then 


(3.4) 


^ [0  1 (n* ) - 61 

*¥r 


- * • W(0,1) 

n'-H» 


where 


’ E{A2}  2E{A1Y1)  E{Y2} 

.e2{a1)  e{a1)e{y1}  e2{y1>_ 

and  6'(n')  is  defined  by  (2.4). 

Proof: 

Define  7^  = Aj  - CYj  and  = Aj  - 0Yj  for  j - 1,2 N(n'), 

and  i = 1,2,  ....  n'  . Then 


(3.5)  M'  = E{NX }M  = E{N1 >62 


(3.6) 


>n'  [0  * (n' ) - Q| 


«N,)  I"'z- 

E{ Y^  } /^M7 


EtY^n' 


t 


where  E{.\  } 1 by  Lemma  2.4.  Also 


(3.7) 


e{n  ) 

1 ^1  l 

E{Y  } /ITm' 


T 

i 


N(n')V(Z1) 


1 

'EO^) 


Ol/2 


zMpH 
N(n')  J 


+ 


n' 

l *1 

K(n')+1  1 

’ V(Zi>l 

n'  

E(N1)  . 

1/2 

wliere  K(n') 


N(n') 

l Nj  and  ViZ^  = E{A2)  - 20E{A1Y1}  + 02E{Y2}  . 


> i 


Since  0 < n'  - K(n')  < • Lew®*  2.2  end  the  SLLN  yield 


(3>8)  ETCT[f(^+'~'ll(n')"'?]  T 1 ' 

x n ***° 


Using  a similar  argument  for  the  positive  and  negative  parts  of  Zj 
and  as  was  used  in  Theorem  2.3, 


(3.9) 


n'  / 

l vjtf  - 

•\4.1  17 


K(n')+1 


0 w.p.l  . 


By  the  central  limit  theorem  for  sums  of  iid  random  variables, 


(3.10) 


N(n') 

zi 


I(n')V{Z^}  n’->® 


W(0,1)  , 


since  is  positive  and  finite  by  the  assumed  conditions  on 

A1  and  Y^  . Combining  (3.7)  - (3.10)  and  applying  Theorem  4.  A.  8 

from  Chung  (1968)  f If  K JL  K*  and  L -*►  1 w.p.l  , then  K L JL  K*1 
L n n n n 

n-x®  n-x°  n-Ko  J 


yields 


(3.11) 


E{N  } l z; 

1 1 1 p 

E{Y^}  /n'H'  n'-*« 


N(0,1)  . 


(3.12) 


EfY^n' 

1 w.p.l 


E(N.)  I Y! 

1 


from  (2.9),  a second  application  of  Chung's  Theorem  4.4.8  to  (3.11) 
and  (3.12)  yields  the  desired  result. 

Note  that  (3.2)  and  (3.5)  are  the  same  when  U ■ {u } and 
V - (v)  . Also  note  that  Theorem  2.5  Implies  Theorem  2.3. 

2.4  The  Ratio  Kstimators 

Let  ((X^.Y^)}^  be  n iid  observations  of  the  bivariate  random 

2 

variable  (X,Y),  which  has  means  yx  and  viy  , variances  and 

2 

Oy  , and  covariance  o^y  . Assume  that  these  moments  are  finite  and 
that  Y,  > 0 w.p.l.  Estimator*  of  0 • yv/p  are  called  ratio 

estimator a. 

The  simple  ratio  estimator  of  0 (already  encountered  in 
(1.6)]  is 

(4.1)  0c(n)  - , 

Y(n) 

n 

where  X(n)  - £ X./n  and  Y(n)  * ][.  Y,/n  are  the  sample  means. 

1 11 

For  sufficiently  large  n , a power  series  expansion  shows  that  the 
dominant  bias  of  (4.1)  is  of  order  1/n  , 

(4.2)  E{0t.  (n)}  3 0 + a/n  + b/n^  + ... 

where  a and  b are  constants  of  the  expression.  As  a result,  an 
extensive  literature  has  developed  in  an  effort  t(.  create  efficient, 
estimators  of  0 which  reduce  the  bias  of  (4.2). 

One  such  estimator  is  the  jackknife  ratio  estimator  due  originally 
to  Quenouille  (1944,  1956)  and  Durbin  (1959)  [see  Miller  (1974)  for 
a comprehensive  review  of  the  hi*  tory  and  properties  of  this  estimator]. 
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The  essence  of  the  jackknife  is  to  divide  the  n observations 

A 

into  g groups  of  h observations  each  (n  ■ gh)  . Let  0£(n) 

A 

be  the  estimate  of  6 based  on  the  full  sample,  and  define  6_^(g,h) 
to  be  the  corresponding  estimator  based  on  a sample  of  size  (g  - l)h  , 
where  the  ith  group  of  size  h has  been  deleted.  The  generalized, 
jackknife  estimator  of  0 is 


(4.3) 

where 


0 (g.h)  - £ l Mg.h)  , 
J 8 i«l  1 


(4.4)  61(g,h)  - g0c(n)  - (g  - l)0_i(&.!  ) , i - 1,2 


are  called  pseudo-values  by  Tukey  (1958). 

Since  (4.3)  is  not  restricted  to  ratio  estimation,  it  is  inter- 
esting to  consider  the  problem  of  estimating  the  mean  of  a sequence  of 

iid  random  variables  Z,,L Z . Then 

1 . n 


ih 


(4.5) 


0j (g»h)  * £ l Z f i - 1,2 

(i-l)h+l  k 


so  that  (4.3)  is  a generalization  of  a batch  means  estimator. 

When  originally  proposed,  Quenouille  considered  the  "simple" 
jackknife  with  g * 2 and  found  that  the  technique  eliminated 
the  1/n  term  from  any  bias.  Since  the  (X  ,Y  ) are  iid,  the 
®i(g.h)  are  interchangeable  random  variables.  Hence,  it  is  straight- 
forward to  show  from  (4.2)  that 


(4.6) 


E{0j(B,h)}  = 0 - b/[g(g  - l)ln  ] + 


■ -•  - — . n...  — — ■ — 
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2 

which  shows  chat  che  bias  co  order  1/n  Is  minimized  for  g ■ n 

2 

(hence,  h ■ 1).  The  1/n  (and  higher)  bias  terms  may  be  eliminated 
by  (repeatedly)  jackknifing  the  pseudo-values  (4.4)  [Schucany,  Gray 
and  Owen  (1971) ) . 

The  principle  attraction  of  the  jackknife  estimator  is  the  ease 

A 

with  which  estimates  of  the  variance  of  0j  can  be  made.  In  an 
abstract,  Tukcy  (1958)  proposed  that  the  pseudo-values  could  be  treated 
as  g approximately  iid  observations  from  which  one  could  expect 


(4.7) 


to  be  approximately  distributed  as  a Student  - t random  variable  with 
g - 1 degrees  of  freedom,  and,  in  the  general  setting  of  U - 
statistics,  Arvesen  (1969)  proved  this  conjecture  for  the  limiting 
case  as  h ♦ « . In  the  same  paper,  Arvesen  also  proved  that  (4.7), 
with  g = n(h  = 1)  , has  a limiting  unit  Normal  distribution  as 
n -►  00  . 


As  a simple  example  of  Arvesen' s theorems,  consider  the  real- 
valued function  f defined  on  (X,Y)  such  that  E{f(X,Y)}  - 6 , 

and  let  {(X^,Y^)}”  be  a sequence  of  iid  observations.  From  his 

2 „ .2 

x < 00  and  oy 


2 2 

Theorem  8,  if  g = n(h  » 1)  , o„  < 00  and  ay  < <*>  , then 


(4.8) 


</g  (Oj  - e) 


8-1 


i=l 


<°i  - V 


172 


g-KO 


/o,e2 

r 2 
°X 

2aXY 

+ il\ 

2 " 

\ 

LPX 

X Y 

“v  J/ 

-Aid 


Also,  under  the  same  conditions,  a special  case  of  his  Theorem  9 yields 
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(4.9) 


2 P 

g-K» 


2°XY 

Vr 


9 


where  5 denotes  convergence  in  probability.  Thus,  0 has  the 

J 

A 

same  asymptotic  variance  as  e^n)  (see  (3.1)  in  section  2.3).  The  form 

of  (4.9)  suggests  the  jackknife  variance  estimator  of  0 , 

J 


(4.10) 


v(eJ(g,h)}  - 


g(g  - 1) 


i=l 


<°i  - 


°J>‘ 


The  alternative  to  using  the  simple  or  jackknife  estimators  of 
0 is  to  use  a ratio  estimator  such  as  Beale's  ratio  estimator 
[Beale  (1962)], 


(4.11) 


0B(n) 


Ksi 

Y(n) 


’l  + -1- 


XY 


n X(n)Y(n) 


1 + - 
■ n 


YY 


Y2(n) 


or  Tin's  ratio  estimator  [Tin  (1965)) 


(4.12) 


9T(n) 


X (n) 
Y(n) 


1 + - I 


XY 


YY 


n\x(n)Y(n)  Y2(n) 


where 


(4.13) 


‘XY  " n - 1 £ ^Xi  " X (n)  ] [ Y - Y(n)] 

i**l 


is  the  sample  covariance.  The  corresponding  estimator  of  the 
variance  of  (4.11)  and  (4.12)  [also  (4.1)]  is  the  siirple  variance 
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estimator 


V{0(n)}  - -jj 

X(n)  ' 

2 

8 XX 

28xy 

-YY 

Y(n) 

X2(n) 

X(n)Y(n) 

*■  -2 
YZ(n) 

Since  there  is  a positive  probability  that  (4.14)  can  be  negative 
for  small  samples,  the  jackknife  ratio  estimator  and  its  associated 
variance  estimator  are  usually  preferred  over  (4.1),  (4.11),  (4.12), 
and  (4.14).  This  choice  is  further  justified  by  the  growing  emperical 
evidence  (Rao  and  Webster  (1966),  Hutchison  (1971),  and  Iglehart 

A 

(1974)]  which  suggests  that  the  mean-square-error  of  0j  for  small 

A 

samples  is  smaller  than  that  of  6 and  is  about  the  same  as  those 

c 

A A 

of  0B  and  0T  . 

The  simple  ratio  estimator  (4.1)  and  the  jackknife  ratio  estimator 
(4.3),  along  with  their  associated  variance  estimators  (4.14)  and  (4.10), 
will  be  used  for  the  emperical  results  of  Chapters  3 and  4.  Since  the 
almost  regenerative  method  can  provide  more  observations  than  the 
regenerative  method,  for  a fixed  run  length,  one  might  expect  that 
the  estimators  (4.1),  (4,3),  (4.10),  and  (4.14)  would  perform  better 
with  the  almost  regenerative  data  than  with  the  regenerative  data. 

These  issues  will  be  discussed  in  the  next  two  chapters. 


i 
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CHAPTER  3 

APPLICATION  OF  THE  ALMOST  REGENERATIVE  METHOD 
TO  A TANDEM  QUEUEING  NETWORK 


In  this  chapter  the  potential  power  of  the  almost  regenerative 
method  (ARM)  will  be  demonstrated  by  simulating  a queueing  network 
consisting  of  two  M/M/1  queues  in  tandem.  The  ARM  is  compared  with 
the  regenerative  method  (RM),  the  fixed  time  increment  method  (FTM), 
imd  a third  method  called  the  random  observer  athod  (ROM) . This 
last  method  consists  of  taking  observations  of  the  process  each  time 
a Poisson  "observer"  arrives.  The  measure  of  comparison  is  the  mean- 
square-error  (MSE)  of  the  estimator  of  total  expected  delay  in  queue. 
For  the  small  sample  size  considered,  the  ARM  estimates  of  delay 
are  20. A,  0.6,  and  23.7  percent  more  efficient  (in  terns  of  MSE) 
than  the  estimates  generated  by  the  RM,  FTM,  and  ROM,  respectively, 
when  using  the  simple  ratio  estimator.  When  using  the  jackknife 
ratio  estimator,  the  estimated  improvements  in  efficiency  are  42.5, 
3.0,  and  25.4  percent,  respectively.  In  addition,  the  estimates  of 
the  variance  of  the  estimator  of  delay  are  superior  for  the  ARM 
generated  data. 

Similar  emperical  results  are  shown  for  simulations  of  an 
H4/H4/I  tandem  queue  network,  where  represents  a hyperexponential 
distribution  with  squared  coefficient  of  variation  of  4. 


3. 1 The  Almost  Regenerative  State  Space  Approach 


Consider  a system  which  consists  of  two  M/M/1  queues  in  tandem 


where  the  departure  stream  of  the  first  queue  is  the  arrival  stream 
for  the  second  queue.  Assume  that  arrivals  to  the  first  queue  occur  . 


1 

I 


s 


t 


K 


at  Poisson  rate  X and  service  at  each  queue  Is  exponential  at  rato 
U . Assume  further  that  each  queue  is  stable  in  the  sense  that 
p = X/w  < 1 . Let  be  the  steady  state  delay  (not  including 
service)  of  an  arbitrary  customer  in  queue  k , k ■ 1,  2,  and  let 
his  total  delay  in  the  network  be  + D2  • Then  it  is 

well-known  that 


(1.1) 


d = E{Dt)  - E{D1>  + E{D2)  , 


where 


(1.2) 


E{V 


ae2{ sk> 

1 - p 


k 


1,2 


and  EtS^)  “ 1/y  is  the  expected  service  time  in  system  k . 

Suppose  an  arrival  finds  the  network  in  state  (i,J)  , where  1 
is  the  number  of  customers  in  system  1 and  j is  the  number  In 
system  2.  From  a workload  point  of  view,  a job  is  merely  a package 
of  "work"  for  the  server.  When  i,j  >0  , the  total  work,  V , in 
the  network  found  by  an  arrival  consists  of  2i  - 1 service  times 
for  those  customers  in  the  first  system,  j - 1 additional  service 
times  for  those  customers  in  the  second  system,  plus  the  remaining 
service  times  of  the  two  customers  in  service.  From  the  "memoryless" 
property  of  the  exponential  distribution,  the  expected  remaining 
service  time  of  a customer  in  service  is  the  same  as  the  expected 
service  time  of  an  arbitrary  customer.  Thus,  the  expected  work 
in  the  network  found  by  an  arrival  who  finds  state  (i,j)  is 
E{V  | (i , j) ) = (2i  + j)E{S}  . 

Since  service  and  interarrival  times  in  the  M/M/1  tandem  network 


i 
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are  exponentially  distributed,  it  is  well-known  that  the  epochs  at 
which  arrivals  find  state  (i,J)  , for  fixed  1 and  j , form  a 
regenerative  process.  Let  p^  be  the  steady  state  probability 
that  the  network  is  in  state  (i,J)  . Since  Poisson  arrivals  "see" 
time  averages  [cf.  Wolff  (1970)],  p^  is  also  the  proportion  of 
time  arrivals  find  the  network  in  state  (i,j)  . Hence,  the  arrival 
rate  of  those  arrivals  who  find  state  (i,j)  is  Xp^  , and  the  mean 
time  between  such  arrivals  is 

(i.  3)  • 

lj  APij 

It  is  also  well-known  [Reich  (1957)]  that,  for  the  M/M/1  tandem 
network. 


(1.4) 


(1  - P) 


2 i+j 
P J 


where  p = X/p  e (0,1)  , so  that  (1.3)  has  a unique  minimum  at 
i " j ■ 0 . Thus 


(1.5) 


E{T  } 
oo 


1 

Ml  - p)2 


is  the  smallest  expected  length  of  a regenerative  cycle,  over  all 
possible  choices  of  (i,j)  , for  the  M/M/1  tandem  network.  Since 
E{Toq}  increases  and  is  unbounded  as  p increases  toward  1 , the 
regenerative  method  might  require  excessively  long  simulation  runs 
to  obtain  a specified  degree  of  accuracy. 

The  almost  regenerative  method,  on  the  other  hand,  cun  lessen 
the  severity  of  this  problem.  Define  the  set  of  states 
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(1.6)  A(n)  - {(i,j)  : 2i  + j - n ; i,j  - 0,1,  ...  ) 

for  n - 0,1,  ....  Then  each  time  an  arrival  to  the  tandem  network 
finds  any  state  (i,j)  e A(n)  , for  fixed  n , he  finds  n equivalent 
services  ahead  of  him.  The  expected  work  in  the  network  found  by  an 
arrival  to  the  set  A(n)  is  thus 

I 

(1.7)  E{V  | A(n) } - nE{S)  . 


Arrivals  to  the  set  A(n  ) , for  fixed  n , form  an  "almost"  regener- 

o o 

ative  process.  Note  that  this  process  is  the  superposition  of  a 
number  of  regenerative  processes  [namely  those  generated  by  the 
elements  of  A(nQ)]  , and  thus  the  asymptotic  results  of  Theorems 
2.2  and  2. A will  hold. 

Analogous  to  the  regenerative  case,  define  T'  to  be  the  time 

n 

between  successive  arrivals  who  find  n equivalent  services  in  the 
netwci.i*  ^r.d  to  be  the  steady  state  probability  that  there  are  n 

QO 

equivalent  services  in  the  network.  The  distribution  {p^}g 
then  given  by 

(1*8)  p'  - l (1  - p)2pi+^,  n =»  0,1,  ... 

" A(n) 

which  simplifies  to 

(1. 9)p2n  ■ (1  “ P)pn(l  - Pn+1),  P2n+i  “ PP2n»  n “ 0*1*  •••• 

As  in  (1.3),  the  expected  time  between  successive  arrivals  finding  n 
equivalent  services  in  the  network  is 


■t 


t 


I 

i 


L. 


ifearSv  i-**.  ***** 


32 


i 


(1.10) 


Since  one  purpose  of  the  ABM  Is  to  Increase  the  number  of 

observations,  for  fixed  run  length,  over  that  possible  from  the  RM, 

It  Is  useful  to  know  the  conditions  under  which  (1.10)  Is  less  than 

(1.3).  Straightforward  algebra  reveals  that  (1.10)  achieves  its 

minimum  value  at  n , where  n is  the  nearest  non-negative  integer 

o o 

to 


(1.11) 


ln2p 

lnp 


for  p e (0,1)  . From  (1.11)  nQ  is  positive  if,  and  only  if,  p 

2/3 

is  greater  than  (1/2)  , or  about  .630.  As  a consequence, 


E{T^  ) < E{Too } if,  and  only  if,  the  same  condition  is  satisfied, 
o 

Otherwise  E{T'  } » E{T  } . For  this  reason  a value  for  p 
n oo 

o 

which  is  greater  than  .630  was  used  for  the  simulation  experiment 


described  below. 


3.2  Simulation  of  the  M/M/1  Tandem  Network 


* Ml  qp 
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fixed  time  increment  method,  FTM,  where  the  conitant  sampling  interval 

size  is  A ■ E{T'  } . The  fourth  method  consists  of  taking  observations 
n 

o 

at  random  Intervals  according  to  a Poisson  process  with  mean  E{T'  } . 

no 

The  purpose  of  this  "random  observer"  method,  denoted  by  ROM,  is  tc 

demonstrate  the  obvious  fact  that  bservations  of  the  cumulative 

delay  process  may  not  be  independent,  even  when  the  times  between 

observations  are.  In  contrast,  however,  observations  of  the  cumulative 

process  for  the  ARM  [A(n  )]  will  be  shown  to  act  as  if  they  were 

o 

iid. 


The  system  utilization  factor, 
queue,  with  X = 1 and  p - 1.25. 

(2.1)  E{Dt)  - 

And  from  (1.3),  (1.10)  and  (1.11), 

i - j - 0 
o Jo 

(2.2) 


p , was  set  at  .8  for  each  M/M/1 
From  (1.1)  and  (1.2) 

6.400  . 


, E{T  } - 25  , 
oo 


n *=  4 , E{T'  } - 16  . 
o n 

o 

One  hundred  independent  replications,  of  m sample  observations 
each,  were  made  for  each  of  the  four  methods:  the  RM(0,0),  the 

ARM[A(4)  ] , the  FTM(A  = 16),  and  the  ROM  with  mean  inter-observation 
time  E{T^}  = 16  . The  parameter  m was  set  equal  to  100  for  each 
method,  with  the  exception  of  the  RM(0,0)  for  which  m was  set 
according  to 

(2.3)  m - 100  E{T'  }/E{T  } , 

n oo 

o 

so  that  m ® 64  for  the  101(0,0)  • The  expected  number  of  arrivals 
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(160,000)  to  the  syetem  ves  then  the  seme  for  each  method.  The 
Interarrival  and  service  times  vers  generated  using  the  CDC  6400 
RAUF  pseudo-random  number  generator.  The  initial  conditions  for 
each  replication  were  the  same;  l.e.,  no  customers  in  the  system 
and  all  servers  Idle. 

The  following  notation  will  be  helpful  in  understanding  the 
simulation  results.  For  each  simulation  run,  r replication  esti- 

^ j*  * 

mates  of  d , {d^  , and  r replication  estimates  of  V(dj}, 

A * j» 

(V(dj)}^  , were  collected.  The  estimators  for  the  jth  replication 
were  based  on  m observations  of  the  cumulative  delay  process, 

{A,. i » and  m observations  of  the  cumulative  arrival  process, 
^Nij^i»l  ^see  Law  (1974) J.  The  estimates  of  interest  are  the  estimate 
of  d , 

(2.4)  i - i J 'i 

j-1  3 


the  estimated  variance  of 


» 


(2.5) 


*2 

°d(J) 


the  estimate  of  the  standard  deviation  of  d , 


(2.6) 


<7  ’ 


the  estimated  MSE  of  d^  , 


MSE(d 


,)  ■ 7 I W,  - d): 

J 1 j-1  J 


(2.7) 
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and  the  estimated  standard  deviation  of  MSE(dj), 


(2.8) 


MSE 


r 11/2 

7(TTTT  {\Mi  - d)2  ' MSEW))2 

j-1  ■» 


The  two  estimates  of  d are  the  simple  ratio  estimator  and  the 
jackknife  ratio  estimator  (see  Chapter  2,  section  2.4).  The  estimated 
variance  of  is  an  external  estimate  since  it  is  formed  at  the 

end  of  the  simulation  run.  An  internal  estimate  of  the  variance  of 
dj  may  be  obtained  from  the  m observations  { (A^  ,1^)  usin8 

the  simple  variance  estimator  and  the  jackknife  variance  estimator 
described  in  Chapter  2,  section  2.4.  These  internal  estimates  will 

*\t 

be  referred  to  with  tildes  (^) , od(j)  » with  corresponding  estimated 
standard  deviation  given  by 


(2.9) 


'V 


V£ 

/r 


The  simulation  results  are  summarized  in  Tables  3.1  and  3,2. 

When  using  the  simple  ratio  estimator,  the  ARM  estimates  of  the 
MSE(d)  are  20.4,  0.6,  and  23.7  percent  less  than  the  estimated 
MSE(d)  generated  by  the  RM,  FTM,  and  ROM,  respectively.  These 
improvements  in  efficiency  are  42.5,  3.0,  and  25.4  percent,  respect- 
ively, when  using  the  jackknife  ratio  estimator.  The  standard  deviation 
of  the  ARM  estimates  of  d and  the  MSE  are  generally  the  smallest 
among  tiie  four  methods,  thus  lending  additional  confidence  to  the 
results. 

The  jackknife  ratio  estimator  appears  to  be  less  biased  than 
the  simple  ratio  estimator,  as  expected.  Also,  the  jackknife  estimates 
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of  d from  Che  RM  and  ARM  data  are  nearly  unbiased,  while  the 
jackknife  estimates  of  d from  the  FTM  and  ROM  data  show  no  bias 
improvement  over  the  simple  ratio  estimates  of  the  same  data.  This 
Is  attributable  to  the  fact  that  the  FTM  and  ROM  observations  are  not 
lid,  and  thus  the  bias  of  each  observation  is  not  the  same;  i.e., 
the  bias  for  the  initial  observations  is  less  than  the  bias  for  the 
final  observations.  Note  that  the  jackknife  ratio  estimator  is 
designed  for  equal  bias  in  all  segments  of  the  data.  On  the  other 


hand,  the  RM  observations  are  iid,  and  the  AR11  observations  "act" 
as  if  they  were  ltd.  This  property  of  the  ARM  data  will  be  discussed 


in  the  next  section. 

A 

Assuming  that  the  external  estimates  of  the  variance  of  d^ 
arc  unbiased  (this  is  not  unreasonable  since  this  estimate  is  based 
on  r independent  replications  of  the  experiment),  it  is  interesting 

A 

to  note  from  Table  3.2  that  the  internal  estimates  of  V(d^)  » both 
for  the  simple  variance  estimator  and  for  the  jackknife  variance 
estimator,  are  biased  considerably  more  for  the  RM  data  than  they 
are  for  the  ARM  data.  In  fact,  the  bias  reduction  is  about  80  percent 

when  using  the  ARM  data.  The  fact  that  the  ARM  generates  more  obser- 
vations, which  act  as  if  they  were  iid,  than  the  RM,  for  a fixed 
run  length,  is  primarily  responsible  for  this  bias  reduction.  The 
FTM  and  ROM  internal  estimates  of  V(dj)  are  not  shown  since  the 
FTM  and  ROM  generated  data  is  highly  correlated  and  produces;  variance 
estimates  of  poor  quality. 

To  give  additional  confidence  to  the  M/M/1  tandem  network 


results,  the  same  experiment  was  also  conducted  for  the  H,/H./l  tandem 

4 4 

network.  The  symbol  represents  a hyperexponential  distribution 
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function  with  squared  coefficient  of  variation  k . The  squared 

2 

coefficient  of  variation,  CV  (X)  , for  the  random  variable  X , 
with  non-zero  mean,  is  CV^(X)  ■ V(X)/E^(X)  . The  distribution 
for  the  experiment  was  taken  to  be  a mixture  of  the  unit  step 
function  at  zero  and  an  exponential  distribution;  l.e., 

(2.10)  P{X  < x)  = H^x)  - bUQ(x)  + (1  - b)e'6x 

where  U (x)  ■ 1 if  x ■ 0 , and  U (x)  ■ 0 if  x > 0 . When 
o o 

E(X)  and  k are  specified,  b - (k  - l)/(k  + 1)  and  0 ■ 

2[(k  + l)E(X)]  * . In  the  case  of  an  arrival  process,  this  distri- 
bution represents  batch  Poisson  arrivals,  where  the  batches  occur 
at  rate  0 and  the  size  of  the  batch  is  geometrically  distributed. 
In  the  case  of  service,  represents  batch  service,  where  the 

service  rate  for  a b;.tch  is  0 . 

As  in  the  M/M/1  tandem  network  simulation,  p was  set  at  .8 
for  each  queue.  Since  the  theoretical  steady  state  distribution 
of  the  total  number  of  customers  in  the  network  is  not  known,  it 
was  estimated  by  making  a pilot  run  of  1500  arrivals  using  the  FTM. 
The  mode  of  the  estimated  distribution,  after  smoothing,  was  near 
8,  and  thus  the  almost  regenerative  data  collection  set  [defined  by 
(1.6)]  was  chosen  to  be  A(8). 

The  results  for  the  H^/H^/l  tandem  model  are  summarized  and 
compared  in  Tables  3.3  and  3.4.  To  facilitate  comparisons,  the 
assumption  was  made  that  the  regenerative  estimates  were  unbiased. 
This  assumption  is  actually  not  true  since  the  jackknife  estimators 
are  biased  even  for  iid  observations.  However,  it  provides  a good 
approximation  since  the  bias  contribution  is  small  in  the  MSE 


estimates  (see  Table  3.3).  As  in  the  M/M/1  case,  the  ARM  dominates 

A 

the  other  methods  for  estimating  d and  V(d^)  . 

3.3  Testa  of  the  Renewal  Hypothesis 

The  tests  of  Appendix  1 were  applied  to  the  almost  regenerative 
cycle  length  process  {Y^}*'  , where  is  the  time  between  the 
(i  - l)st  and  the  ith  arrivals  to  find  the  system  in  state  A(nQ)  , 
and  the  corresponding  cumulative  delay  process  {Aj}!j|  » where  A^ 
is  the  cumulative  network  delay  during  Y^  , for  the  M/M/1  and  the 
H^/H^/l  tandem  queue  models.  The  purpose  of  the  tests  .as  to  deter- 
mine whether  the  observed  processes  act  like  those  of  a regenerative 
process. 

Figure  3.1  contains  a plot  of  the  estimated  autocorrelation 
functions,  based  on  1000  observations,  for  the  M/M/1  cycle  length 
process  and  the  M/M/1  cumulative  delay  process  generated  by  the 
ARM[A(4)].  With  the  exception  of  the  spike  at  lag  14,  there  appears 
to  be  no  significant  autoregressive  behavior  in  either  process. 

The  correlation  at  lag  14  is  difficult  to  explain;  possibly  it  is 
due  to  slight  nonrandom  behavior  in  the  CDC  RANF  pseudo-random 
number  generator  which  is  known  to  have  short  periods  in  its  low 
order  bits  [Hutchinson  (1966)].  Tests  of  RANF,  using  the  same 
initial  seeds  as  for  the  simulation  experiment,  confirmed  that  some 
periodicities  do  exist  in  its  estimated  spectrum.  It  is  interesting 
to  note  that  the  estimated  autocorrelation  functions  for  the  RM(0,0), 
applied  to  the  M/M/1  tandem  network,  also  showed  what  appeared  to 
be  significant  autoregressive  behavior  at  lags  9 and  34.  The  un- 
explained autoregressive  behavior  persisted,  but  at  different  lags, 
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AUTOCORRELATION  FUNCTION  FOR  25  LAGS 
M/M/1  TANDEM  ARM  WITH  A ( N ) STATE  SPACE 

FIGURE  3.1 
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when  new  random  number  seeds  vers  used  for  RAMF.  It  appears  that 
this  spurious  correlation  Is  not  significant  for  the  purposes  of  the 
present  Investigation. 

The  estimated  autocorrelation  function  for  the  M/M/1  FTM  cumulative 
process,  shown  in  Figure  3.2,  has  significant  correlation  for  at 
least  five  lags.  The  corresponding  processes  for  the  H^/H^/l  tandem 
model,  shown  in  Figures  3.3  and  3.4,  lead  to  similar  conclusions. 

Table  3.5  contains  a summary  of  the  renewal  test  statistics  for 
each  of  the  observed  processes.  It  is  interesting  that  the  processes 
generated  by  the  ARM  comfortably  pass  all  tests,  at  the  95  per  cent 
significance  level,  for  both  the  M/M/1  and  the  H^/II^/1  tandem  models. 
Further,  none  of  these  processes  exhibits  any  discernable  periodicitcs, 
at  the  95  per  cent  significance  level,  in  its  estimated  spectrum  (not 
shown).  Thus,  the  ARM,  when  applied  to  the  M/M/1  and  the  H^/H^/l 
tandem  network  models,  generates  observations  which  act  as  though 
they  were  iid.  As  expected,  the  FTM  processes  do  not  pass  any  of  the 
renewal  tests  at  the  95  per  cent  level. 


3.4  Distribution  of  the  Times  Between  Transitions 

In  this  section  expressions  will  be  derived  for  the  times  between 
state  transitions,  conditioned  on  the  number  of  equivalent  services 
in  the  system,  for  the  M/M/1  tandem  network,  and  it  will  be  shown 
that  these  times  are  distributed  "almost"  exponentially  for  certain 
values  of  p ■ X/p  . 

First,  consider  a non-negative  random  variable  X ~ F(t)  ; i.e.. 


i 

I 

* 


s 

« 

'i 

i 

f, 


P{X  < t}  “ F(t)  . The  Laplace-Stieltjcs  transform  (LS-transform)  of 
X is  defined  by 
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SUMMARY  OF  TEST  STATISTICS  FOR  THE  TANDEM  QUEUE  MODEL 

TABLE  3.5 


(4.1) 


%(a)  -/  e"f 


'dF(t) 


for  s > 0 . The  LS-transform  can  be  Interpreted  as  an  expectatior. 
in  the  following  sense 


(4.2) 


}(s)  - E{e"sX)  . 


This  definition  will  be  useful  in  the  following  discussion. 

Let  T be  the  time  between  state  transitions  in  the  M/ll/1  tandem 
network  model.  In  section  3.1  the  steady  state  probability  there 
are  n equivalent  services  in  the  network  was  found  to  be 

(4.3)  pin  - (1  - p)p"(l  - pn+1),  P2n+1  = Pp2n*  n = °*1 

Since  the  LS-transform  of  an  exponential  distribution  with  parameter 
Y is  y/(y  + s)  . the  conditional  distribution  of  T , given  the  number 
of  equivalent  services  in  the  network,  has  the  LS-transform 


(4.4) 


E{e”sT| 2n} 


X + s 


— + H a + L±J±_ 

X + y + s 2n  X + 2u  <•  s 


a2n) » n m 1 » 2 , ... 


(4.5) 


E{e”sT| 2n  + 1} 


X + y , 

X + y + s 2n+l 


+ L±_2 u__  ,,  . . 

X + 2y  + s (1  b2n+l 
n ■ 0,1,  . . . 


where 
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(4.6) 


a2n  ■ <Po.2»  + pn,o>  [(iiJ  JA(2ll)plj]"1  • n ' 2-2-  - 
b2n+l  ' Po,2n+l  [(liJ)cA{2n+1)Plj]  1 " ’ °>1’  — 


and  ptj  is  the  asymptotic  probability  that  there  are  i customers 
in  system  1 and  j customers  in  system  2.  From  (4.4)  and  (4.5) 
it  is  easy  to  see  that,  w.p.l, 


(4.7) 

and 


E{e“sT|2n}  - (1  - p) 


* + H + 1 + 2u 

A+p+s  P A + 2p  + s 


as  n -*•  “ 


(4.8) 


E{e_sT j 2n  + 1} 


-*■  * + 2u 

A + 2 + s 


as  n 


Thus,  (T | 2n}  converges  to  a hyperexponential  distribution  while 
(T  | 2n  + 1}  converges  to  an  exponential  distribution  as  n becomes 
large.  And,  as  p - 1 and  n - - , the  steady  state  distribution 
of  T is  an  exponential  distribution. 

The  squared  coefficient  of  variation  of  T given  n is  defined  by 


(4.9) 


CV2{t| n}  - JOlhl  m lililn) 
E2{T|n}  E2 { T j n) 


1 » n - 0,1 


From  (4.7)  and  (4.8) 


(4.10) 

and 


CV2(T|2n)  - l + f (1  - p)  , as  „ 


(4.11) 


CV  { T | 2n  + 1}  ■+  1 , as  n » 


More  interesting,  however,  is  the  behavior  of  T , given  n , for 
small  n . From  (4.4)  and  (4.5) 


2(l-p  ) 1 N ■ r/  ±£Z>LLA±IiL  HR  l£TF2— LL~.E U _ j 

[(l-p)(2+p)(l+pn)  + p(l+p)(l-pn“L)]2 


n 


1.2, 


• • • 


and 

(4.13)  CV2(T| 2n+l}  - 2(1  - p11*1) 


[ (l-p)pn(2+p)2  + (l-pn)  (l+p)2J  1 

[ (l-p)pn(2+p)  + (l-pn)(l+p)]2 


» 


n * 0 f X i • • i • 


Table  3.6  shows  selected  values  of  (4.12)  and  (4.13)  as  functions 
of  p and  n . For  all  values  of  p and  n , {Tj2n}  behaves 
very  much  like  an  exponential  random  variable  (which  has  a coefficient 
of  variation  of  1).  However,  {T|2n  + 1}  has  a relatively  high 
coefficient  of  variation  for  small  n > 0 . 

Unconditioning  on  n in  (4.4)  and  (4.5)  yields 


(4.14) 


} 


I+s  + 2p<1"u) 


A+n  2 \+2y 

A+y+s  L A+2y+s 


which  has  squared  coefficient  of  variation 


(4.15) 


CV2{T}  - 5p4  + 2p3  - 19c2  + 12a  + 4 
(P2  - P + 2) 2 


Table  3.6  contains  values  of 


cv2{t} 


as  a function  of  p , from  which 


one  can  argue  that  T is  distributed  "almost"  exponentially  for  large 
values  of  p . It  is  interesting  to  note,  as  shown  earlier,  that  the 
almost  regenerative  method  is  most  powerful  for  values  of  p > (l/2)2^3 


OOOOOOOOO 


S»>r^uS-#^N«SHH 
OOOOOOOOO 


ffs'»OO‘^<*r*'N00 
O »\  rs.  10  cn  «n  .h  «-t 
HOOOOOOOO 


N .H 

H 

sf 

(N 

H 

Q 

o 

n o\ 
H O 
• • 

lA 

O 

• 

(N 

o 

• 

o 

• 

8 

0 

8 

• 

8 

• 

8 

• 

s 

O 

• 

8 

• 

O (N  m H r-  <N  ■-( 

I58SS88 


O rs  m H O 

8^318 


o o o © 

So  o o 

o o o 


or'-p'.»-ta\vO<sio\i/>fNLn 

§^005\MtO(ONNN'} 

OOOOOOOOOO 


rH  I tn 
H 


«Nno\r^r^^»Noaoo 

OOOOOOOOOO 


o\ 

CO 

st 

m 

u-i 

iTl 

m 

5 oo 

o 

H 

CM 

rs 

CN 

CN 

<N 

CM 

©oors.mstu’tmu'immm 

gr^o\oooooooo 

OOiHH^HiHdr-ttHH 


H 

m 

•A 

m 

m 

*» 

8 

• 

8 

• 

*t 

*» 

>t 

St 

st 

*» 

«t 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

OH<sn^ms«r^oo0\o  8 


49 


CHAPTER  4 

USING  THE  ALMOST  REGENERATIVE  METHOD  TO  OBTAIN  ORDER  OF 
MAGNITUDE  IMPROVEMENTS  IN  MSE  EFFICIENCY 


A simple  queueing  network  model  is  presented  in  this  chapter  for 
which  the  ARM  can  provide  a 60  percent  improvement  in  MSE  efficiency 
over  the  RM,  and  a 50  percent  improvement  over  the  FTM,  for  estimating 
the  expected  total  delay  in  queue.  Equally  important,  the  ARM 
generated  data  provides  better  estimates  of  estimator  variance  than 
the  RM.  These  emperical  findings  have  strong  implications  for  the 
potential  power  of  the  ARM  when  applied  to  general  queueing  networks, 
since  it  appears  that  the  relative  efficiency  of  the  ARM  increases 
with  system  complexity. 

In  addition,  comparisons  of  the  simple  ratio  estimator  with  the 
jackknife  ratio  estimator  and  the  simple  variance  estimator  with 
the  jackknife  variance  estimator  show  that  the  jackknife  estimators 
are  preferable  when  making  estimates  from  ARM  generated  data. 


4,1  The  M/M/c  Tandem  Queue  Network 


Let  pn  be  the  steady  state  probability  that  there  are  n 
customers  in  a single  M/M/c  queue,  where  arrivals  occur  at  Poisson 
rate  A and  the  service  time  at  each  channel  is  exponentially  dis- 
tributed at  rate  u . If  p £ A/cu  < 1 , it  is  well-known  that 


(1.1) 

where 


I 


P (A/p)  / n ! , n = 0,1,  ...,  c-i 


| p (A/cp)n  ° , n = c,  c+1,  ... 


i-l 


c 

I 

n”0 


ppQ (A/u) 

(c  - 1 j ! (cii  - / ) 


(1.2) 


n! 


I- 
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It  is  also  well-known  that  the  expected  delay  in  queue  of  a typical 
customer,  given  delay  ia  positive,  is  distributed  exponentially  with 
parameter  (cp  - X)  . Hence 


(1.3) 


E{D) 


Ip.- 


mp0(X/w) 


(c  - l)l(cp  - X) 


2 ‘ 


Now  consider  a network  consisting  of  two  M/M/c  queues  in  tandem, 
so  that  departures  from  the  first  system  are  arrivals  at  the  second 
system.  Again,  assume  that  arrivals  to  the  first  system  are  Poisson 
at  rate  X and  that  the  service  time  at  each  channel,  for  each  queue, 
is  exponentially  distributed  at  rate  p . Let  p^  represent  the 
steady  state  probability  that  there  are  1 customers  in  the  first 
system  and  j in  the  second  system;  i.e.,  the  network  is  in  state 
(i,j).  From  Reich  (1963),  the  departure  stream  from  the  first 
system  is  Poisson  at  rate  X , and  p^  ■ p^  . Also,  the  expected 
total  delay  in  system  is  d = E{D^}  » 2E{D}  . 

One  possible  choice  of  regenerative  state  transition  for  this 
tandem  network  occurs  vheu  an  arrival  to  the  first  system  finds  the 
network  in  state  (i,j),  for  fixed  i and  j . Using  the  same 
arguments  as  in  Chapter  3,  it  is  well-known  that  the  mean  time  between 
such  regeneration  points  is 


(1.4) 


E{V 


Xp 


ij 


This  expression  is  minimized  for  the  state  (i  ,1  ) that  satisfies 

o o 

iQ  “ jQ  and  iQ  < \/p  < i + 1 . As  in  Chapter  3,  the  regenerative 

method  with  regenerative  state  (i  ,1  ) will  be  referred  to  as  the 

o o 
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Consider  the  set  of  states 


(1.5)  A(n1,n2)  - {(i,j)  : i + j c [n^n^  , 1,J  - 0,1,  ...) 

for  n^,  n2  ■ 0,1 An  almost  regeneration  point  will  be  said 

to  occur  whenever  an  arrival  to  the  first  system  finds  any  one  of  the 
states  (i,j)  e A(n^,n2)  , for  fixed  n^,n2  . The  probability  of 
such  an  event  is 


(1.6) 


>'  “ l 

nl,n2  A(ni,n2) 


ij 


U0 

2 n 

I l p 

n=n-  i=0 


i.n  -i  ’ 


and  the  expected  time  between  these  events  is 


(1.7) 


E{T' 

* ^2 


} 


i 


Xp' 

Vn2 


The  almost  regenerative  metliod  with  almost  regenerative  state  A(n^,n2) 
will  be  referred  to  as  the  AKM(A(n^ ,n2) ] . 

4.2  Simulation  of  the  M/M/2  Tandem  Network  Using  ARM  [A(n^,n0)] 

As  shown  emperically  in  Chapter  3,  the  ARM  can  generate  more 
observations,  which  act  as  if  they  were  iid,  than  the  RM  for  the  same 
simulation  run  length.  As  a result,  the  ARM  can  be  more  efficient, 
in  terms  of  the  MSE  of  the  estimator,  than  the  RM  or  the  FTM  when 
estimating  the  expected  delay,  d , for  simple  queueing  networks.  One 
might  expect  that  the  MSE  efficiency  of  the  ARM  could  be  further 
improved  by  increasing  the  number  of  states  in  A(nj,n2)  . But 
increasing  the  size  of  A(n^,n2)  involves  a balance  between  increased 
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observation  frequency  end  more  pronounced  serial  correlation  in  the 

observations.  In  addition,  computer  CPU  time  Increases  as  A(n^,n^) 

becomes  larger  since  more  observations  are  processed  for  the  same 

run  length.  In  this  section  the  ARM(A(n^  tn^)  ] will  be  compared 

with  the  RM(i,1)  and  the  FTM(A-E{T'  })  for  various  values  of 

nl’n2 

n^  and  n^  • 

For  c ■ 2 , A ■ 1 , and  u ■ .625  (so  that  p ■ .8),  (1.4) 

is  minimized  when  i ■ j * 1 . From  (1.1)  - (1.7)  it  is  straight- 

o o 

forward  to  compute 

d - 5.689,  E{T  } - 81.0  , E{T„}  - 31.641, 

oo  ll 

E(T!.}  - 12.361,  E<T!,  1 - 6.180,  E{T.',)  - 4.175, 

(2.1)  44  45  46 

E{T^}  - 3.171,  E{T^y}  - 2.578,  E{T^g)  - 2.202, 

E{T^g}  - 1.933,  E{TJ9)  - 1.733. 


One  hundred  independent  simulation  replications,  k m observations 
each,  were  made  for  each  data  collection  method.  The  parameter  m 
was  set  to  ensure  the  same  expected  number  of  arrivals  to  the  system 
(123,000)  for  each  method.  The  interarrival  and  service  times  (for 
each  system)  were  generated  with  separate  random  number  seeds,  which 
were  the  same  for  each  method,  using  the  CDC  6400  RANF  pseudo-random 
number  generator.  The  estimates  of  d were  based  on  the  cumulative 
delay  process,  {Aj , and  the  cumulative  arrival  process,  {Nj . 

The  parameters  of  interest  were  computed  using  ralations  (2.4)  - 
(2.8)  of  Chapter  3,  section  3.2.  Each  replication  began  with  identical 
initial  conditions  (i.e.,  no  customers  in  queue  and  all  servers  idle). 

The  simulation  results  are  sumnurized  in  Tables  4.1,  4.2,  4.3,  and 
4.4.  Referring  to  Tables  4.1  and  4.2,  the  simple  ratio  estimates  of 
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MSE(dj)  for  Che  ARM  are  up  to  60  percent  less  than  the  RM  simple 


ratio  estimates  of  MSE(d^)  and  up  to  53  percent  less  than  the 


corresponding  FTM  estimates.  This  improvement  reaches  a peak  for 


the  ARM[A(3,8)J.  The  ARM  jackknife  estimates  of  MSE(dj)  are  nearly 


45  percent  less  than  the  RM  jackknife  estimates  of  MSE(d^)  and  up 
to  43  percent  less  than  the  corresponding  FTM  estimates.  This  im- 
provement also  reaches  a maximum  for  the  ARM[A(3,8)].  The  standard 


deviation  of  the  ARM  estimates  of  MSE(dj)  is  generally  less  than 


the  standard  deviation  of  the  UM  and  FTM  estimates,  thus  increasing 
confidence  in  these  results. 

As  observed  for  the  M/M/1  tandem  network  of  Chapter  3,  the 
ARM  and  RM  jackknife  estimates  of  d for  the  M/M/2  network  have 
smaller  biases  than  the  corresponding  simple  ratio  estimates  of  d . 
However,  the  FTM  jackknife  estimates  of  d have  about  the  same 
bias  as  the  corresponding  simple  ratio  estimates.  This  result  is 
directly  attributable  to  the  superior  iid  qualities  of  the  ARM 
generated  data  when  compared  to  the  FTM  data. 


The  internal  and  external  estimates  of  V (d^ ) are  shown  in 


Table  4.3  for  the  R11  and  ARM  generated  data.  The  variance  estimates 
using  the  FTM  data  are  not  shown  since  the  serial  correlation 
between  the  FTM  observations  degrades  the  estimator  accuracy  and 
increases  the  CPU  time  necessary  to  form  the  estimator  (since  Jagged 
autocorrelations  must  be  included  in  the  computation). 


From  Table  4.3,  the  internal  ARM  estimates  of  V(d^)  are  less 


biased  and  move  stable  than  the  RM  estimates,  where  the  measure  of 
stability  is  the  estimated  standard  deviation  of  the  estimator.  It 
is  interesting  to  observe  that  the  RM(0,0)  internal  variance  statistic 
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underestimates  significantly,  but  the  corresponding  standard  deviation 
estimate  does  not  provide  an  indication  of  this  fact;  i.e.,  the 
RM(0,0)  internal  variance  estimate  has  a large  bias  while  its 
standard  deviation  implies  that  the  variance  estimate  is  relatively 
good.  This  problem  is  also  characteristic  of  the  ARM  variance 
estimates,  but  the  severity  of  the  problem  is  reduced.  The  impor- 
tance of  this  improvement  cannot  be  overly  emphasized,  since  in 
most  simulations  the  variance  of  the  estimator  of  interest  is  usually 
an  internal  estimate  based  on  a single  simulation  run. 

When  comparing  the  simple  variance  estimator  and  the  jackknife 
variance  estimator  for  the  different  data  collection  methods  (see 
Table  4.3),  the  internal  jackknife  variance  estimator  is  consistently 
within  one  standard  deviation  of  the  external  estimate.  This  is  not 
always  the  case  for  the  simple  variance  estimator. 

The  renewal  test  statistics  for  the  ARU[A(n^ .n^) ) are  shown  in 
Table  4.4.  The  lack  of  any  discernable  periodicities  in  the  associated 
spectrums  (not  shown),  at  the  95  percent  level,  supports  the  conclusion 
from  Table  4.4  that  the  jackknifed  transformations  of  the  ARM[A(4,4)] 
and  the  ARM[A(4,5)]  generated  observations  act  as  if  they  were  iid. 

The  fact  that  data  generated  by  larger  Afa^^)  sets  do  not  pass 
the  renewal  tests  does  not  seem  to  degrade  the  corresponding  estimates 
of  d and  V(dj)  . It  seems  clear,  however,  that  the  jackknife 
transformation  improves  the  iid  properties  of  the  data. 

There  are  two  reasons  why  the  ARM  estimates  of  d are  more 
efficient  than  the  FTM  estimates  of  d . First,  the  AFuM  generated 
data  has  better  iid  properties  than  the  FTM  data,  and,  secondly,  the 
jackknife  transformation  of  the  data  further  improves  its  iid 
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qualities.  Also,  It  appears  that  the  jackknife  estimators  do  not 
reduce  estimator  bias  when  the  data  has  significant  serial  correlation. 

Increasing  the  size  of  ACn^^)  does  not  appear  to  cause 
severe  degredatlon  of  these  results  when  n^  and  n£  are  within  a 
"reasonable"  range  of  each  other.  However,  as  the  separation  between 
n^  and  n^  increases,  the  number  of  observations,  for  a fixed  run 
length,  also  increases,  thereby  increasing  the  CPU  time  to  process 
the  additional  observations. 

4.3  Comments  on  the  Implied  Power  of  the  Almost  Regenerative  Method 

When  comparing  the  emperical  results  of  this  chapter  with  those 
of  Chapter  3,  it  appeals  that  the  MSE  efficiency  of  the  ARM  estimators, 
relative  to  the  RM  and  the  FTM,  increases  with  system  complexity. 

In  addition,  the  internal  estimates  of  estimator  variance  are  improved 
when  the  estimates  are  based  on  the  ARM  generated  data  rather  than 
the  RM  generated  data.  These  observations  have  far-reaching  conse- 
quences, since  it  is  not  difficult  to  imagine  queueing  netv.'orks  for 
which  the  ARM  could  provide  an  order  of  magnitude  improvement  in 
MSE  efficiency. 

These  properties  of  the  ARM  are  due  primarily  to  the  fact  that 
the  ARM  can  generate  observations  which  act  as  if  they  were  iid, 
unlike  the  FTM,  and  that  the  ARM  can  generate  more  observations,  for 
a fixed  run  length,  than  the  RM.  To  illustrate  this  second  point, 
the  theoretical  expected  tine  between  some  M/M/c  tandem  network 
regeneration  and  almost  regeneration  points  are  shown  in  Table  4.5 
for  selected  values  of  X , u , and  c . 
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p 

c 

RM 

E{T  }(1)  E{T.  , }(2)  E{ T * } 

oo  i i n ,n 

0 0 0 0 

ARM 

(3)  Eli'  )W) 

Vno+1 

.8 

2 

81.0 

31.6 

12.4 

6.2 

.8 

3 

316.8 

38.2 

12.5 

6.2 

.8 

4 

1345.5 

45.0 

12.7 

6.3 

.8 

5 

5929.0 

52.1 

13.5 

6.7 

.9 

2 

361.0 

111.4 

25.6 

13.0 

.9 

3 

1612.0 

121.3 

25.9 

3 3.0 

.9 

4 

7892.5 

130.5 

26.0 

13.0 

.9 

5 

40671.6 

139.3 

26.2 

13.1 

.95 

2 

1521.0 

421.3 

53.0 

26.5 

.95 

3 

7237.7 

438.8 

53.1 

26.6 

.95 

4 

37996.4 

454.3 

53.1 

26.6 

.95 

5 

210786.8 

468.5 

53.2 

26.7 

(1) 


Expected  time  between  arrivals  who  find  state  (i,j)  * (0,0)  . 


Expected  time  between  arrivals  who  find  state  (i  ,i  ),  where 

i ' j and  i < X/y  < i + 1 . 
o o o o 

(3) 

Expected  time  between  arrivals  who  find  some  state  (i,j)  e A(n  ,n  ), 

o o 

where  n minimizes  E{ T * } . 

o n ,n 

(4)  o ° 

Expected  time  between  arrivals  who  find  some  state  (i,1)  e A(n  ,n  xl) 

o o 

where  n minimizes  E{T'  . 

o n ,n  +1 

o o 
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CHAPTER  5 

THE  ALMOST  REGENERATIVE  METHOD  IN  PERSPECTIVE 

There  are  two  related,  but  distinct,  Issues  Involved  when 
conducting  a simulation  experiment.  The  internal  problem  is  concerned 
with  the  conditions  under  which  data  is  collected  from  the  simulation, 
and  the  external  problem  is  concerned  with  the  use  made  of  the  data 
to  form  reliable  estimates  of  the  parameters  of  interest. 

The  regenerative  method  (RM)  solves  the  internal  problem  for 
simulations  of  regenerative  processes  by  generating  observations 
which  are  iid.  As  a result,  the  external  problem  is  eased  since 
autocovariances  need  not  be  computed  when  estimating  the  variance 
of  the  estimators.  However,  the  11M  does  have  at  least  one  serious 
failing.  It  cannot,  be  applied  to  simulations  of  stochastic  processes 
which  are  not  regenerative,  nor  can  it  be  applied  effectively  to 
regenerative  processes  which  have  excessively  long  expected  times 
between  regeneration  times.  A second,  less  serious,  defect  concerns 
the  external  problem  of  forming  ratio  estimators  of  the  data.  All 
known  ratio  estimators  are  biased,  and  the  bias  depends  on  the 
regenerative  state  chosen. 

The  alternative  methods  for  collecting  simulation  data  arc  the 
standard  fixed  time  increment  method  (FTI1),  which  makes  no  use  of 
any  prior  knowledge  of  the  process  under  consideration,  and  approximate 
regeneration  methods.  The  approximate  methods  can  be  classified  as 
cither  state  space  transformation  methods,  which  distort  the  original 
process  by  inducing  a regenerative  approximation  to  the  process  of 
interest,  and  the  almost  regenerative  methods. 


i 


The  purpose  of  this  dissertation  hus  been  to  develop  An  almost 
regenerative  method  which  alleviates  the  internal  data  collection 
problem  created  by  the  RM  (i.e.,  possibly  infrequent  observations) 
and  yet  exploits  the  desireable  property  of  the  RM  for  producing 
lid  observations. 

When  compared  to  t lie  RM,  the  ARM  has  been  shown  emperically  to 
have  the  potential  for  an  order  of  magnitude  reduction  in  the  MSE 
of  the  total  delay  in  queue  estimator  for  simple  queueing  networks. 

In  fact,  this  relative  improvement  increases  with  system  complexity. 
Similar  results  were  shown  when  comparing  the  ARM  with  the  FTM. 
Moreover,  the  ARM  variance  estimates  have  a smaller  bias  and  a 
smaller  estimates  standard  deviation  than  the  RM  variance  estimators. 
The  programming  considerations  necessary  to  implement  the  ARM  are 
no  more  complicated  than  those  required  to  implement  the  PJI  or  the 
FTN  for  most  event-oriented  simulations. 

However,  the  ARM  is  not  witiiout  fault.  The  choice  of  the  almost 
regenerative  pair  of  sets  (U,V)  , which  trigger  data  collection, 
remains  largely  an  art,  since  an  intelligent  choice  of  (U,V) 
depends  on  the  characteristics  of  the  stochastic  process  being 
simulated,  the  parameters  being  measured,  and  the  desired  estimator 
accuracy.  Increases  in  the  sizes  of  U and  V provide  more  obser- 
vations for  fixed  run  length,  but  only  at  the  possible  expense  of 
variance  estimator  accuracy.  That  is,  once  the  data  loses  its 
similarity  to  iid  data,  the  autocorrelation  between  observations 
must  be  included  in  the  variance  estimates  and  this  requires  an 
additional  pass  through  the  data,  and  probabily  a degradation  of  the 


1 


variance  estimator. 
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As  an  external  issue,  the  Jackknife  transformation  seems  to 
improve  the  "lid-like"  properties  of  the  data  (see  Table  4.3, 

Chapter  4),  thus  improving  variance  estimates.  This  fact  needs 
to  be  balanced  with  the  larger  storage  requirements  of  the  jackknife 
estimator,  when  compared  to  the  storage  requirements  of  the  other 
estima  tors. 


■i 


\ 
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APPENDIX  I 


TESTING  FOR  RENEWAL  PROCESSES 


The  tests  in  this  Appendix  can  be  found  in  Cox  and  Lewis  (1966). 

Let  {X^}^  be  a random  sample  for  a stationary  random  variable 

X with  unknown  continuous  distribution  Fv(t)  , and  let  F (t)  be 

X o 

a completely  specified  distribution.  The  tests  of  this  Appendix  arc 
designed  to  test  the  null  hypothesis 


li  : Fv  - F 

o X o 


against  the  two-sided  alternative 


Hi:  * * • 

i X o 


The  first  two  tests  which  will  be  used  later  are  distribution-free 
in  the  sense  that  the  null  hypothesis  does  not  depend  on  the  specific 
form  of  Fv  and  F . This  property  is  a consequence  of  the  pro- 
bability  integral  transformation  which  states  that  if  the  random 
variable  X has  the  distribution  F , then  U * F(X)  is  uniformly 
distributed  over  [0,1]  . Thus,  the  probability  integral  transformation 
reduces  the  problem  to  one  of  testing  whether  the  observations 
{FQ(Xi)}^  are  uniformly  distributed. 

The  first  test  is  the  well-known  Komogorov-Smirnov  test.  Define 
the  statistic 


where 


Cn  - sup  | Fn<t)  - FQ(t)|  - max  (Dn,D~) 

— 00<  t <°° 

D*  - sup  (F  (t)  - F (t) } » max  F (X,.,)l 

n „ ..T  n o ..  ( n o (i)J 

1<“  i«l, . . . ,n 


v*v  asmok  Wflir  **(^* 


Dn  - sup  {F  (t)  - F 

-»<t«»  ° 


In  these  expressions,  Fn(t)  Is  the  eraperical  distribution  obtained  by 


number  of  X's  < t 

F (t)  - L_^L_ 

n n 


and  X(1)  is  the  ith  order  statistic  for  the  sequence  {X^J  . If 
dn>a  Is  a value  such  that 


P{D  > d } - a , 
n n,a  ’ 

then  the  null  hypothesis  is  rejected  at  level  a if  the  observed  value 
of  is  greater  than  d^  . It  is  well-known  [cf.  Darling  (1957)) 

that  Dn  and  Dn  have  the  asymptotic  distribution 


(A.  1) 


lira  P{D+/n  > d}  - e~2d 
n-x»  n 


for  d c [0,«)  ; (A. 1)  is  generally  valid  for  n > 20  . 

The  second  type  of  test  Is  based  on  the  Cramdr-von  Mises  statistic 


- Fo(t)J2dFo(t) 


which  has  the  computational  form 


Again  the  null  hypothesis  is  rejected  at  level  a If  the  observed  value 


* 1 
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of  W2 
n 


is  greater  than  w where 
n,a 


w } 
n,a 


a . 


Both  of  these  tests  are  measures  of  the  distance  between  the  two 
distributions.  It  is  interesting  to  note  that  these  two  tests  are 
distribution  free,  consistant,  and  do  not  require  arbitrary  groupings 
of  the  data.  Other  better  known  tests,  for  example  the  chi-squared 
goodness  of  fit  test,  do  not  have  all  these  properties.  Many  other 
less  powerful  tests  for  uniformity  are  in  Knuth  (1969). 


A.l  Serial  Correlation  Coefficient  Tests 


There  are  essentially  two  types  of  tests  for  renewal  processes: 
those  based  on  the  serial  correlation  coefficients,  and  those  based 
on  the  spectrum  of  intervals.  The  serial  correlation  coefficient 
tests  will  be  described  first. 

As  before,  let  {X^}"  be  a sequence  of  observations  X^  + X 

00  00 

from  some  stationary  stochastic  process  {X^q  where  X . The 

usual  unbiased  estimator  for  the  serial  correlation  coefficient  of 
lag  j » » is  given  by 


(A.  3) 


where 


R 


j 


i n-j 

n E IXi  " x<n>HXi+j  - X(n)J 


j » 0,1,  . . . , n-1 


is  an  estimate  of  the  covariance  at  lag  j , 

Cj  - E((Xi  - EX][Xi+j  - EX])  , 


■* r> ."v *, 


for  all  i,j  - 0,1,  ....  Since 


r,  - r n 

--  ■ — * W<0,1)  w.p.l  . 

/iT^T 


[Bartlett  (1955)],  a test  at  the  a level  that  r^  - 0 is  to  reject 


independence  if 


Q 

lril  ’7 


'n  - 1 


where  Qo/2  is  the  upper  a/2  point  of  the  unit  Normal  distribution. 
Although  this  test  is  simple,  it  has  a number  of  shortcomings . First, 
it  is  difficult  to  extend  this  test  to  sample  serial  correlation 
coefficients  of  lag  greater  than  one.  Secondly,  the  small  sample 
distribution  of  r^  under  the  independence  hypothesis  is  not  known; 
clearly,  it  depends  on  the  distribution  of  times  between  events. 

A test  designed  to  alleviate  this  second  difficulty  is  the  Wald- 
Wolfowitz  rank  product-moment  test  for  lag  1 [Wald  and  Wolfowitz  (1943)] 
The  central  idea  is  to  replace  the  X^'s  In  (A. 3)  by  their  ranks, 
where  the  rank  of  is  the  argument  of  its  corresponding  order 

statistic.  Evidently,  this  procedure  has  the  advantage  of  controlling 
the  Influence  of  spurious  observations.  Wald  and  Wolfowitz  showed  that 
the  resulting  rank  product-moment  statistic,  r^  , is  asymptotically 
normally  distributed  with  mean  and  variance  given  by 


E{V  " TI  (n  * 1)(n  + 1)(3n  + 2>  • 


V{r  ) - 5°6  ± 16n5  - 14nA  - SOn3  - 35n2  + 64n  + 44 


720(n  - 1) 
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A.  2 Tests  Based  on  the  Spectrum  of  Inf  rvals 

More  complicated!  but  more  powerful,  tests  for  renewal  processes 
are  based  In  the  spectrum  of  Intervals.  Consider  the  stationary 

oo 

sequence  {X^^  which  has  a spectral  density  defined  by 


(A.  4) 


f(u) 


1_ 

2s 


l r.  coa(ju) 

ja-a  J 


-71  < U < TT 


and  a power  spectrum  defined  by 


2 1 

(A. 5)  0 f(w)  » -j-  l C COs(ju>)  , -7T  < u < 71 

J— • J 

2 

where  o is  the  variance  of  the.  r.v.  X . If  the  {X^}  are  iid,  then, 
from  (A. 4), 

(A. 6)  f(w)  ■ , -77  < w < Tt 


so  that  a test  for  renewal  process  is  to  compute  the  spectral  density 
function  and  check  to  see  if  it  is  of  the  form  (A.  6).  Define 
u>p  * 2irp/n  , where  p ■ 1,2,  nj  . It  can  be  shown  that  the 

periodogramt  defined  by 


I (u  ) ■ -r — 
n p 2im 


J-l 


V 


-iu  j 
P 


a)  i4  0 (mod  2ti) 
P 


where  i «*  *^1  , is  an  estimator  of  the  power  spectrum  (cf.  Brillinger 
(1975)].  Bartlett  (1955)  and  Hannan  (1960)  showed,  among  other  things, 
that  when  the  X^  are  iid  the  periodogram  asymptotically  has  an  ex- 
ponential distribution  with  mean 


p*' 


I tftmy’iflff***  f'Wl Wf fV  **$?  1 ^ **'*$*’  . 


(A.  8) 


O2  2 

E^In^Wp^  " 2n  " 0 f^“p)  » **»p  *<  0 (mod  2u)  . 
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Cox  and  Lewis  (1966)  suggest  testing  the  quantities 


(A.  9) 


(1) 


Ci  vy/ty 
Ci  W'**-.) 


as  samples  from  a uniform  distribution  to  test  the  hypothesis  that 
In(u,p)  is  exponentially  distributed.  The  Kolmogorov-Smirov  and 
Crara6r-von  Mises  statistics  are  particularly  well-suited  for  this  kind 
of  test. 


A. 3 A Test  of  the  Tests 

The  following  linear  model  was  simulated,  for  different  values  of 
a and  i , 

(A*  10)  xt  ■ «xt_£  + et  , ct  ~W(0,1)  , t - 0,  +1,  + 2,  .... 

to  test  the  sensitivity  of  the  tests  in  this  section  for  accepting  the 
renewal  hypothesis.  The  results  of  these  tests  are  summarized  in 
Table  A.l.  In  all  those  cases  where  the  renewal  hypothesis  was 
accepted  on  the  basis  of  the  tests,  strong  pcriodicies  appeared  in  the 
estimated  spectrum  of  {Xt)  at  the  95  per  cent  significance  level. 
Thus,  the  renewal  hypothesis  was  rejected  in  each  case. 

Let  X(n)  - l Xx/n  be  the  usual  estimator  of  the  mean  of  the 
process  defined  by  (A. 10).  Since  this  process  is  stationary,  one 
may  easily  show  that  the  variance  of  X(n)  is  given  by  [of.  Fishman 
(1973b),  pg.  280) 


■ 

J 


tot 


TESTS  OF  THE  RENEWAL  PROCESS  TEST  STATISTICS 


I 

! 


If  one  were  to  assume  that  no  serial  dependence  existed  in  n 
successive  observations  of  the  process  defined  by  (A.  10),  then  the 

percent  of  the  true  variance  given  by  (A.  11)  that  would  be  lost 
can  be  expressed  as 


l 

< 

* 
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